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Abstract 

We define a class of "algebraic" random matrices. These are random matrices for which the Stieltjes 
transform of the limiting eigenvalue distribution function is algebraic, i.e., it satisfies a (bivariate) poly- 
nomial equation. The Wigner and Wishart matrices whose limiting eigenvalue distributions are given by 
the semi-circle law and the Marcenko-Pastur law are special cases. 

Algebraicity of a random matrix sequence is shown to act as a certificate of the computability of 
the limiting eigenvalue density function. The limiting moments of algebraic random matrix sequences, 
when they exist, are shown to satisfy a finite depth linear recursion so that they may often be efficiently 
enumerated in closed form. 

In this article, we develop the mathematics of the polynomial method which allows us to describe 
the class of algebraic matrices by its generators and map the constructive approach we employ when 
proving algebraicity into a software implementation that is available for download in the form of the 
RMTool random matrix "calculator" package. Our characterization of the closure of algebraic probability 
distributions under free additive and multiplicative convolution operations allows us to simultaneously 
establish a framework for computational (non-commutative) "free probability" theory. We hope that the 
tools developed allow researchers to finally harness the power of the infinite random matrix theory. 

Key words Random matrices, stochastic eigen- analysis, free probability, algebraic functions, resul- 
tants, D-finite series. 



1. Introduction 

We propose a powerful method that allows us to calculate the limiting eigenvalue distribution of a 
large class of random matrices. We see this method as allowing us to expand our reach beyond the well 
known special random matrices whose limiting eigenvalue distributions have the semi-circle density [38] , the 
Marcenko-Pastur density [18], the McKay density [19] or their close cousins [8,25]. In particular, we encode 
transforms of the limiting eigenvalue distribution function as solutions of bivariate polynomial equations. 
Then canonical operations on the random matrices become operations on the bivariate polynomials. We 
illustrate this with a simple example. Suppose we take the Wigner matrix, sampled in Matlab as: 

G = sign(randn(N))/sqrt(N) ; A = (G+G' ) / sqrt (2) ; 

whose eigenvalues in the N — > oo limit follow the semicircle law, and the Wishart matrix which may be 
sampled in Matlab as: 

G = randn(N,2*N)/sqrt(2*N) ; B = G*G'; 
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whose eigenvalues in the limit follow the Marccnko-Pastur law. The associated limiting eigenvalue distribu- 
tion functions have Stieltjes transforms m^jz) and TOb(z) that are solutions of the equations L^ z (m, z) = 
and L^ z (m, z) = 0, respectively, where 

L^ z (m, z) = m 2 + zm + 1, L^(m, z) = m 2 z - (-2 z + 1) m + 2. 

The sum and product of independent samples of these random matrices have limiting eigenvalue distribution 
functions whose Stieltjes transform is a solution of the bivariate polynomial equations L^ B (m, z) = and 
L AB (m,z) = 0, respectively, which can be calculated from L^ z and L^ z alone. To obtain L^t B (to, z) we 
apply the transformation labelled as "Add Atomic Wishart" in Table[7]with c = 2, p\ = 1 and Ai = 1/c = 0.5 
to obtain the operational law 

L2+ B (m,z) = Lt ^ m , z -—L—j . (1.1) 

Substituting L^L = m 2 + z m + 1 in and clearing the denominator, yields the bivariate polynomial 

L^+ B (m,z) = m 3 + {z + 2) m 2 - (-2 2+ l)m + 2. (1.2) 

Similarly, to obtain we apply the transformation labelled as "Multiply Wishart" in Table[7]with c = 0.5 
to obtain the operational law 

L^f (to, z) = Lt ( (0.5 - 0.5zm) to, Z ) . (1.3) 
\ 0.5 — 0.5zm / 

Substituting L^ z = m 2 + zm+ 1 in (jl.3p and clearing the denominator, yields the bivariate polynomial 

L^ B (m,z) = m 4 z 2 — 2 m 3 z + m 2 + A mz + 4. (1.4) 

Figure [T] plots the density function associated with the limiting eigenvalue distribution for the Wigner and 
Wishart matrices as well as their sum and product extracted directly from L^ B (m, z) and L^ B (m, z). In 
these examples, algebraically extracting the roots of these polynomials using the cubic or quartic formulas 
is of little use except to determine the limiting density function. As we shall demonstrate in Section [STJ 
the algcbraicity of the limiting distribution (in the sense made precise next) is what allows us to readily 
enumerate the moments efficiently directly from the polynomials L^ B (to, z) and L^f (m, z). 



1.1 Algebraic random matrices: Definition and Utility 

A central object in the study of large random matrices is the empirical distribution function which is 
defined, for an N x TV matrix An with real eigenvalues, as 

A Number of eigenvalues of Ajv < x 



For a large class of random matrices, the empirical distribution function F An (x) converges, for every x, 
almost surely (or in probability) as TV — > oo to a non-random distribution function F A (x). The dominant 
theme of this paper is that "algebraic" random matrices form an important subclass of analytically tractable 
random matrices and can be effectively studied using combinatorial and analytical techniques that we bring 
into sharper focus in this paper. 
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(a) The limiting eigenvalue density function for the GOE and Wishart ma- 
trices. 




(b) The limiting eigenvalue density function for the sum and product of 
independent GOE and Wishart matrices. 

Figure 1: A representative computation using the random matrix calculator. 
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Figure 2: A random matrix calculator where a sequence of deterministic and stochastic operations performed 
on an algebraic random matrix sequence A at produces an algebraic random matrix sequence Bat. 
The limiting eigenvalue density and moments of a algebraic random matrix can be computed 
numerically, with the latter often in closed form. 



Definition 1 (Algebraic random matrices). Let F A {x) denote the limiting eigenvalue distribution function 
of a sequence of random matrices A.n- If a bivariate polynomial L mz (m, z) exists such that 

m A {z)= I —^—dF A (x) zeC+\R 



is a solution of L raz {mA{z) 1 z) = then A^v is said to be an algebraic random matrix. The density function 
fA '■= dF A (in the distributional sense) is referred to as an algebraic density and we say that Ajy £ -M a ig, 
the class of algebraic random matrices and fA € V a ig, the class of algebraic distributions. 



The utility of this, admittedly technical, definition comes from the fact that we are able to concretely 
specify the generators of this class. We illustrate this with a simple example. Let G be an n x m random 
matrix with i.i.d. standard normal entries with variance l/m. The matrix W(c) = GG' is the Wishart 
matrix parameterized by c = n/m. Let A be an arbitrary algebraic random matrix independent of W(c). 
Figure [2] identifies deterministic and stochastic operations that can be performed on A so that the resulting 
matrix is algebraic as well. The calculator analogy is apt because once we start with an algebraic random 
matrix, if we keep pushing away at the buttons we still get an algebraic random matrix whose limiting 
eigenvalue distribution is concretely computable using the algorithms developed in Section HO 

The algebraicity definition is important because everything we want to know about the limiting eigenvalue 
distribution of A is encoded in the bivariate polynomial L^ z (m, z). In this paper, we establish the algebraicity 
of each of the transformations in Figure [2] using the "hard" approach that we label as the polynomial 
method whereby we explicitly determine the operational law for the polynomial transformation L^ z (m, z) 1— > 
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L^ z (m, z) corresponding to the random matrix transformation A i— > B. This is in contrast to the "soft" 
approach taken in a recent paper by Anderson and Zeitouni [3, Section 6] where the algebraicity of Stieltjes 
transforms under hypotheses frequently fulfilled in RMT is proven using dimension theory for nocthcrian 
local rings. The catalogue of admissible transformations, the corresponding "hard" operational law and 
their software realization is found in Section [67J This then allows us to calculate the eigenvalue distribution 
functions of a large class of algebraic random matrices that are generated from other algebraic random 
matrices. In the simple case involving Wigner and Wishart matrices considered earlier, the transformed 
polynomials were obtained by hand calculation. Along with the theory of algebraic random matrices we also 
develop a software realization that maps the entire catalog of transformations (see Tables [7] -0 into symbolic 
Matlab code. Thus, for the same example, the sequence of commands: 

>> syms m z 

>> LmzA = m~2+z*m+l; 

» LmzB = m~2-(-2*z+l)*m+2; 

>> LmzApB = AplusB (LmzA, LmzB) ; 

>> LmzAtB = At imesB( LmzA, LmzB ) ; 

could also have been used to obtain L^ B and L^f. We note that the commands AplusB and AtimesB 
implicitly use the free convolution machinery (see Sectioning to perform the said computation. To summarize, 
by defining the class of algebraic random matrices, we are able to extend the reach of infinite random matrix 
theory well beyond the special cases of matrices with Gaussian entries. The key idea is that by encoding 
probability densities as solutions of bivariate polynomial equations, and deriving the correct operational laws 
on this encoding, we can take advantage of powerful symbolic and numerical techniques to compute these 
densities and their associated moments. 



This paper is organized as follows. We introduce various transform representations of the distribution 
function in Section [2U We define algebraic distributions and the various manners in which they can be 
implicitly represented in [37J and describe how they may be algebraically manipulated in |47J The class of 
algebraic random matrices is described in Section [5~J where the theorems are stated and proved by obtaining 
the operational law on the bivariate polynomials summarized in Section [10 Techniques for determining the 
density function of the limiting eigenvalue distribution function and the associated moments are discussed 
in Sections [771 and [871 respectively We discuss the relevance of the polynomial method to computational free 
probability in Section [97J provide some applications in Scction flO.l and conclude with some open problems in 
Section QX] 



We now describe the various ways in which transforms of the empirical distribution function can be 
encoded and manipulated. 

2.1 The Stieltjes transform and some minor variations 

The Stieltjes transform of the distribution function F A (x) is given by 



1.2 Outline 



2. Transform representations 





(2.1) 



The Stieltjes transform may be interpreted as the expectation 



m,A(z) = 



1 



'X 



x — z 
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with respect to the random variable x with distribution function F A (x). Consequently, for any invertiblc 
function h(x) continuous over the support of dF A (x), the Stieltjes transform rriA(z) can also be written in 
terms of the distribution of the random variable y = h(x) as 



ha(z) = E x 



h(-»(y)-z 



(2.2) 



where ^ (■) is the inverse of h(-) with respect to composition i.e. h(h^ ^ (x)) = x. Equivalently. for 
y = h(x), we obtain the relationship 



E,, 



y-z 



E T 



h(x) 



The well-known Stieltjcs-Pcrron inversion formula [1] 



1 



Ja{x) = dF (x) = — lim Im itia(x + it;). 

7T £-+0+ 



(2.3) 



(2.4) 



can be used to recover the probability density function Ja{x) from the Stieltjes transform. Here and for the 
remainder of this thesis, the density function is assumed to be distributional derivative of the distribution 
function. In a portion of the literature on random matrices, the Cauchy transform is defined as 



1 



-dF A {x) forzeC^Y 



9a{z) -- 

The Cauchy transform is related to the Stieltjes transform, as defined in (|2.1[) . by 

9a{z) = ~m A {z). 



(2.5) 



2.2 The moment transform 

When the probability distribution is compactly supported, the Stieltjes transform can also be expressed 
as the series expansion 

1 °° M A 

^) = ---Ej7TT> ( 2 - 6 ) 

i=i 

about z = oo, where M A := J x^dF A (x) is the j-th moment. The ordinary moment generating function, 
Ha{z), is the power series 



3=0 



A z\ 



(2.7) 



1 



with M A = 1. The moment generating function, referred to as the moment transform, is related to the 
Stieltjes transform by 

(2.8) 



(2.9) 



Ha{z) = —m A - 



The Stieltjes transform can be expressed in terms of the moment transform as 



m A {z) = -~Ha 
z 



The eta transform, introduced by Tulino and Verdu in [32], is a minor variation of the moment transform. 
It can be expressed in terms of the Stieltjes transform as 



Va(z) = -^ m A 



(2.10) 
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while the Sticltjcs transform can be expressed in terms of the eta transform as 



m A (z) = ~T)a (-^). (2.11) 

2.3 The R transform 

The R transform is defined in terms of the Cauchy transform as 

r A (z) = g { A 1} (z)--, (2.12) 

z 

where g A ^ (z) is the functional inverse of g A (z) with respect to composition. It will often be more convenient 
to use the expression for the R transform in terms of the Cauchy transform given by 

r A (g)=z(g)--. (2.13) 
9 

The R transform can be written as a power series whose coefficients are known as the "free cumulants." 
For a combinatorial interpretation of free cumulants, see [28]. Thus the R transform is the (ordinary) free 
cumulant generating function 

oo 

r A (g) = Y / Kf+i9 J - (2-14) 

3=0 

2.4 The S transform 

The S transform is relatively more complicated. It is defined as 



1 + z 



*iW = — T^'W (2.15) 



z 



where T A (z) can be written in terms of the Sticltjcs transform m A (z) as 

T A (z) = --m A (l/z)-l. (2.16) 

z 

This definition is quite cumbersome to work with because of the functional inverse in (|2.15[) . It also places 
a technical restriction (to enable series inversion) that ^ 0. We can, however, avoid this by expressing 
the S transform algebraically in terms of the Stieltjes transform as shown next. We first plug in T A (z) into 
the left-hand side of (|2. 15|) to obtain 

s A (T A {z)) = z. 

This can be rewritten in terms of m A (z) using the relationship in (|2. 16[) to obtain 

. 1 . , , . zm(llz) 
8 A (— m(l/z) - 1) - 



z m(l/z) + z 

or, equivalently: 

, v > m(z) 

s A (-zm(z)-l) = j^—. 2.17 

z m{z) + 1 

We now define y(z) in terms of the Stieltjes transform as y(z) = —zm(z) — 1. It is clear that y(z) is an 
invertible function of m(z). The right hand side of (|2.17p can be rewritten in terms of y{z) as 

. . ss m(z) m(z) 

s A {y{z)) = y- = 7 \L—. 2.18 

y(z) zm(z) + l 
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Equation (|2.18p can be rewritten to obtain a simple relationship between the Stieltjes transform and the S 
transform 

m A (z) = -ys A (y). (2.19) 
Noting that y = —zm(z) — 1 and m(z) = —ys A (y) we obtain the relationship 

y = zys A (y) - 1 

or, equivalently 

-^V. (2.20) 

ysA(y) 

3. Algebraic distributions 

Notation 3.1 (Bivariate polynomial). Let L uv denote a bivariate polynomial of degree D u in u and D v in 

v defined as 

D u D v D u 

l uv = l uv (; •) = E E c ^ ul yk = E u3 - 

3=0 k=0 j=0 

The scalar coefficients Cjk are real valued. 

The two letter subscripts for the bivariate polynomial L uv provide us with a convention of which dummy 
variables we will use. We will gencrically use the first letter in the subscript to represent a transform of the 
density with the second letter acting as a mnemonic for the dummy variable associated with the transform. 
By consistently using the same pair of letters to denote the bivariate polynomial that encodes the transform 
and the associated dummy variable, this abuse of notation allows us to readily identify the encoding of the 
distribution that is being manipulated. 

Remark 3.2 (Irrcducibility). Unless otherwise stated it will be understood that L uv (u,v) is "irreducible'' in 
the sense that the conditions: 

• lo(v), . . . , fe u (i | ) have no common factor involving v, 

• disc L (u) ^ 0, 

are satisfied, where disci(u) is the discriminant of L uv (u,v) thought of as a polynomial in v. 
We are particularly focused on the solution "curves," Ui(v), . . . , ujy u (v), i.e., 

L uv (u,v) =Z Du (u) JJ (u - Ui{v)) . 

i=l 

Informally speaking, when we refer to the bivariate polynomial equation L nv (u,v) = with solutions ut(v) 
we are actually considering the equivalence class of rational functions with this set of solution curves. 

Remark 3.3 (Equivalence class). The equivalence class of L uv (u,v) may be characterized as functions of 
the form L uv (u,v)g(v)/h(u,v) where h is relatively prime to L uv (u,v) and g(v) is not identically 0. 
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A few technicalities (such as poles and singular points) that will be catalogued later in Section ^2 remain, 
but this is sufficient for allowing us to introduce rational transformations of the arguments and continue to 
use the language of polynomials. 

Definition 3.4 (Algebraic distributions). Let F(x) be a probability distribution function and f(x) be its 
distributional derivative (here and henceforth). Consider the Stieltjes transform m(z) of the distribution 
function, defined as 

m{z)= I dF(x) forzeC + \R. (3.2) 



If there exists a bivariate polynomial L mz such that L mz (m(z), z) = then we refer to F(x) as algebraic 
(probability) distribution function, f(x) as an algebraic (probability) density function and say the f € V a ig- 
Here Vaig denotes the class of algebraic (probability) distributions. 

Definition 3.5 (Atomic distribution). Let F(x) be a probability distribution function of the form 

K 

F { x ) = ^Pi^lXuoo), 
i=l 

where the K atoms at Xi g R have (non-negative) weights pi subject to ^2 t Pi = 1 and I[a;,oo) * s the indicator 
(or characteristic) function of the set [x, oo). We refer to F(x) as an atomic (probability) distribution 
function. Denoting its distributional derivative by f{x), we say that f{x) £ V a tom- Here V atom denotes the 
class of atomic distributions. 

Example 3.6. An atomic probability distribution, as in Definition \3. 51 has a Stieltjes transform 

K 



K 

Pi 



, Xi - z 

which is the solution of the equation L rnz (m, z) = where 

K K K 



i—1 i—1 j^ki 

i=i 

Hence it is an algebraic distribution; consequently Vatom C V a \g- 
Example 3.7. The Cauchy distribution whose density 

1 



fix) 



n{x 2 + 1)' 

has a Stieltjes transform m(z) which is the solution of the equation L mz (m, z) = where 

L mz (m, z) = (z 2 + l) m 2 +2zm+l. 
Hence it is an algebraic distribution. 

It is often the case that the probability density functions of algebraic distributions, according to our def- 
inition, will also be algebraic functions themselves. We conjecture that this is a necessary but not sufficient 
condition. We show that it is not sufficient by providing the counter-example below. 
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Counter-example 3.8. Consider the quarter- circle distribution with density function 



f(x) = ^LJ^L for xe [0,2]. 

TT 

Its Stieltjes transform : 

4-2V-^ 2 +41nf- 2+ ^+ 7 ") + 
771(2) = 



2ir 

is clearly not an algebraic function. Thus f(x) ^ T- \i g . 

3.1 Implicit representations of algebraic distributions 

We now define six interconnected bivariate polynomials denoted by L mz , L gz , L rg , L sy , L^ z , and L . We 
assume that L nv (u, v) is an irreducible bivariate polynomial of the form in (|3.1j) . The main protagonist of the 
transformations we consider is the bivariate polynomial L mz which implicitly defines the Stieltjes transform 
771(2) via the equation L mz (m, z) = 0. Starting off with this polynomial we can obtain the polynomial L gz 
using the relationship in (|2.5p as 

L &z {g,z) = L mz (-g,z). (3.3) 

Perhaps we should explain our abuse of notation once again, for the sake of clarity. Given any one polynomial, 
all the other polynomials can be obtained. The two letter subscripts not only tell us which of the six 
polynomials we are focusing on, it provides a convention of which dummy variables we will use. The first 
letter in the subscript represents the transform; the second letter is a mnemonic for the variable associated 
with the transform that we use consistently in the software based on this framework. With this notation in 
mind, we can obtain the polynomial L rg from L gz using (|2.13p as 

L Ig (r,g)=L s Jg,r+^j . (3.4) 

Similarly, we can obtain the bivariate polynomial L sy from L mz using the expressions in (|2.19p and (|2.20[) to 
obtain the relationship 

Based on the transforms discussed in SectinO we can derive transformations between additional pairs of 
bivariate polynomials represented by the bidirectional arrows in Figure [3] and listed in the third column of 
Table [3] Specifically, the expressions in (|2.8[) and (|2. 1 1[) can be used to derive the transformations between 
L mz and L^ z and L mz and L rjZ respectively. The fourth column of Table [3] lists the Matlab function, imple- 
mented using its Maple based Symbolic Toolbox, corresponding to the bivariate polynomial transformations 
represented in Figure [3] In the Matlab functions, the function irreducLuv(u, v) listed in Table [1] ensures 
that the resulting bivariate polynomial is irreducible by clearing the denominator and making the resulting 
polynomial square free. 

Example: Consider an atomic probability distribution with 

F(s)=0.5I [0> oo)+0.5I [llOo) , (3.6) 

whose Stieltjes transform 

0.5 0.5 

m{z) = + , 

— 2 1—2 

is the solution of the equation 

m(0 - 2)(1 - z) - 0.5(1 - 22) = 0, 
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Legend 

m(z) = Stieltjes transform 
g(z) = Cauchy transform 
r(g) = R transform 
s(y) = S transform 
H(z) = Moment transform 
r)(z) = Eta transform 



Figure 3: The six interconnected bivariate polynomials; transformations between the polynomials, indicated 
by the labelled arrows, are given in Tabled 



or equivalently, the solution of the equation L mz (m, z) — where 

L mz (m, z) = m(2 z 2 - 2 z) - (1 - 2z). 

We can obtain the bivariate polynomial L gz (g,z) by applying the transformation in 
polynomial L mz given by (|3.7[) so that 



L gz (g,z) = -g(2z 2 -2z)-(l-2z). 
Similarly, by applying the transformation in (|3.4p wc obtain 



r+ -g)- 2 ( r+ -g 



r+~ 



(3.7) 

to the bivariate 
(3.8) 

(3.9) 



which, on clearing the denominator and invoking the equivalence class representation of our polynomials (see 
Remark l3.3|) . gives us the irreducible bivariate polynomial 



L Tg (r 7 g) = -l + 2gr 2 + (2-2g)r. 
By applying the transformation in (13. 5| to the bivariate polynomial L mz , we obtain 



(3.10) 



- 1-2 



V + 1 



sy \ sy J I \ sy 

which on clearing the denominator gives us the irreducible bivariate polynomial 

L^(s,y) = (l + 2y)s-2-2y. 



(3.11) 



Table [2] tabulates the six bivariate polynomial encodings in Figure [3] for the distribution in (|3.6p , the semi- 
circle distribution for Wigner matrices and the Marcenko-Pastur distribution for Wishart matrices. 



The polynomial method 



12 



Procedure 


Matlab Code 


Simplify and clear the denominator 
Make square free 
Simplify 


function Luv = irreducLuv(Luv,u, v) 

L = numden(simplif y (expand(Luv) ) ) ; 
L = Luv / maple ( 'gcd' ,L,diff (L,u) ) ; 
L = simplify (expand(L) ) ; 
L = Luv / maple ( 'gcd' ,L,diff (L,v) ) ; 
Luv = simplify (expand (L) ) ; 



Table 1: Making L uv irreducible. 



(a) The atomic distribution in Il3,6|l . 



L 


Bivariate Polynomials 




m{2z' 2 -2z)-{l-2z) 


L gz 


-g(2z 2 -2z) - {l-2z) 


L rg 


-1 + 2gr 2 + (2-2g)r 


L sy 


{l + 2y)s-2-2y 




(-2 + 2z)n + 2~ z 




(2z + 2)r)~2-z 



(b) The Marfienko-Pastur distribution. 



L 


Bivariate Polynomials 


Lmz 
L S z 

L rg 
L sy 


czm 2 — (1 — c — z)m + 1 

czg 2 + (1 — c — z) g + I 

(cg-l)r + l 

(cy + l)s-l 

H zc — (zc + 1 — z) fj, + 1 

rj 2 zc + {—zc + 1 — z) 7] — 1 



(c) The semi-circle distribution. 



L 


Bivariate polynomials 


L m z 


m 2 + m z + 1 


L e z 


9 2 - 9 z + 1 


L Tg 


r — g 


L sy 


s 2 y-l 


Lfj.z 


H 2 z 2 — (J, + 1 


L V z 


z 2 j] 2 -77 + I 



Table 2: Bivariate polynomial representations of some algebraic distributions. 
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Label 


Conversion 


Transformation 


MATLAB Code 


I 


^mz "^gz 




function Lmz = Lgz2Lmz(Lgz) 
syms mgz 

Lmz = subs(Lgz,g,-m); 


£gz = L nm(~9, 2) 


function Lgz = Lmz2Lgz(Lmz) 
syms mgz 

Lgz = subs(Lmz,m,-g); 


II 


^gz ^rg 


^gz = ^rg(« ~ 


function Lgz = Lrg2Lgz(Lrg) 
syms rgz 

Lgz = subs(Lrg,r,z-l/g); 
Lgz = irreducLuv(Lgz,g,z); 


^rg -^ g zVy= ' * ) 


function Lrg = Lgz2Lrg(Lgz) 
syms rgz 

Lrg = subs(Lgz,g,r+l/gJ; 
Lrg = irreducLuv(Lrg,r,g); 


III 


^mz ^rg 


^mz ^ ^gz ^ ^rg 


function Lmz = Lrg2Lmz(Lrg) 

syms m z r g 

Lgz = Lrg2Lgz(Lrg); 

Lmz = Lgz2Lmz(Lgz); 


function Lrg = Lmz2Lrg(Lmz) 
syms m z r g 
Lgz = Lmz2Lgz(Lmz); 
Lig = Lgz2Lig(Lgz); 


IV 




L m = L *y( zm + i> zm X ) 


function Lmz = Lsy2Lmz(Lsy) 
syms m z s y 

Lmz = subs(Lsy,s,m/(z*m+l)); 
Lmz = subs(Lmz,y,-z*m-l); 


L w -L mz ( ys, sy ) 


function Lsy = Lmz2Lsy(Lmz) 
syms m z s y 
Lsy = subs(Lmz,m,-y*s); 
Lsy = subs(Lsy,z,(y+l)/y/s); 


V 






function Lmz = Lmyuz2Lmz(Lmyuz) 

syms m myu z 

Lmz = subs(Lmyuz,z,l/z); 

Lmz = subs(Lmz,myu,-m*z); 

Lmz = irreducLuv(Lmz,m,z); 




function Lmyuz = Lmz2Lmyuz(Lmz) 
syms m myu z 
Lmyuz = subs(Lmz,z,l/z); 
Lmyuz = subs(Lmyuz,m,-myu*z); 
Lmyuz = irreducLuv(Lmyuz,myu,z); 


VI 




L mz = L ,,J- zm ' --) 


function Lmz = Letaz2Lmz( Letaz) 
syms m eta z 

Lmz = subs(Letaz,z,-l/z); 
Lmz = subs(Lmz,eta,-z*m); 
Lmz = irreducLuv(Lmz,m,z); 




function Letaz = Lmz2Letaz(Lmz) 
syms m eta z 

Letaz = subs(Lmz,z,-l/z); 
Letaz = subs(Letaz,m,z*eta); 
Letaz = irreducLuv(Letaz,eta,z); 



Table 3: Transformations between the different bivariate polynomials. As a guide to Matlab notation, the 
command syms declares a variable to be symbolic while the command subs symbolically substitutes every 
occurrence of the second argument in the first argument with the third argument. Thus, for example, the 
command y=subs(x-a,a,10) will yield the output y=x-10 if we have previously declared x and a to be 
symbolic using the command syms x a. 
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4. Algebraic operations on algebraic functions 

Algebraic functions are closed under addition and multiplication. Hence we can add (or multiply) two al- 
gebraic functions and obtain another algebraic function. We show, using purely matrix theoretic arguments, 
how to obtain the polynomial equation whose solution is the sum (or product) of two algebraic functions 
without ever actually computing the individual functions. In Section 14.21 we interpret this computation 
using the concept of resultants [31] from elimination theory. These tools will feature prominently in Section 
[37] when we encode the transformations of the random matrices as algebraic operations on the appropriate 
form of the bivariatc polynomial that encodes their limiting eigenvalue distributions. 



4.1 Companion matrix based computation 

Definition 4.1 (Companion Matrix). The companion matrix C a M to a monic polynomial 

a(x) = ao + ai x + . . . + a n -i x"^ 1 + x n 

is the n x n square matrix 



C„ 



(x) 



-a 
-ai 



-ai 



.0 1 — On-l. 

with ones on the sub-diagonal and the last column given by the negative coefficients ofa(x). 



Remark 4.2. The eigenvalues of the companion matrix are the solutions of the equation a(x) = 0. This is 
intimately related to the observation that the characteristic polynomial of the companion matrix equals a(x), 
i.e., 

a{x) = det(.Tl„ - C o(a; )). 

Consider the bivariatc polynomial L uv as in (|3.ip . By treating it as a polynomial in u whose coefficients are 
polynomials in v, i.e., by rewriting it as 



L uv (u,v) = i 
j'=o 



(4.1) 



we can create a companion matrix CJJ V whose characteristic polynomial as a function of u is the bivariatc 
polynomial L . The companion matrix CJJ V is the D u x D u matrix in Table 2] 





Matlab code 


■ -lo(v)/ln u (v) - 

1 -h(v)/l Dn (v) 

'■• -h(v)/ln a (v) 
.0 1 -Z Du -i(v)/Zd u («) . 


function Cu = Luv2Cu(Luv,u) 
Du = doublc(maple('degree',Luv,u)); 
LDu = maple('coeff ',Luv,u,Du); 
Cu = sym(zeros(Du))+ .. 

+diag(ones(Du-l,l),-l)); 
for Di = 0:Du-l 

LtuDi = maple('coeff ',Lt,u,Di); 

Cu(Di+l,Du) = -LtuDi/LDu; 

end 



Table 4: The companion matrix CJj v , with respect to u, of the bivariate polynomial L uv given by (14.11) . 
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Remark 4.3. Analogous to the univariate case, the characteristic polynomial of C]J V is det(ul — CJJ V ) = 
L uv (u, v)//d u (^) Du - Since ln n (v) is not identically zero, we say that det(itl — CJJ V ) = L uv (u,v) where the 
equality is understood to be with respect to the equivalence class of L uv as in Remark \3.S\ . The eigenvalues 
o/C|J v are the solutions of the algebraic equation L uv (u,v) = 0; specifically, we obtain the algebraic function 
u(v). 



Definition 4.4 (Kroncckcr product). If A m (with entries a^) is an m X m matrix and B„ is an n x n 

matrix then the Kronecker (or tensor) product of A m and B„ ; denoted by A m ® B„ ; is the mn x ran matrix 
defined as: 

B 

.^ml^Tl ■ • ■ a mn B T 

Lemma 4.5. I/a, and /3j are the eigenvalues of A m and B„ respectively, then 

1. oti + /3j is an eigenvalue of (A m (X) I„) + (Ln <8 B„), 

,2. [3j is an eigenvalue o/A m <X> B„, 
fori = l,...,m,j = l,...,n. 

Proof. The statements are proved in [16, Theorem 4.4.5] and [16, Theorem 4.2.12]. □ 

Proposition 4.6. Let u\{v) be a solution of the algebraic equation L\ v {u,v) = 0, or equivalently an eigen- 
value of the D\ x D\ companion matrix C^. Let U2{v) be a solution of the algebraic equation L^ v (u, v) = 0, 
or equivalently an eigenvalue of the x companion matrix C"^. Then 

1. uj,(v) = ui(v) + U2(v) is an eigenvalue of the matrix CJJJ = (CJJ^, <E> I_og) + (I_di ® CJJ^), 

2. uz(v) = ui(v)u2(v) is an eigenvalue of the matrix CJJiJ, = CJJJ, <X> C"^.. 

Equivalently 1*3(1") is a solution of the algebraic equation = where L^ v = det(ul — C"^). 

Proof. This follows directly from Lemma 14.51 □ We represent the binary addition and multiplication 



operators on the space of algebraic functions by the symbols EB U and E3 U respectively. We define addition 
and multiplication as in Table [5] by applying Proposition 14.61 Note that the subscript 'u' in EB U and K u 
provides us with an indispensable convention of which dummy variable we are using. Table [5] illustrates 
the EB and M operations on a pair of bivariate polynomials and underscores the importance of the symbolic 
software developed. The (D u +1) x (D v +1) matrix T uv lists only the coefficients cy for the term u % w J in 
the polynomial L uv (u,v). Note that the indexing for i and j starts with zero. 

4.2 Resultants based computation 

Addition (and multiplication) of algebraic functions produces another algebraic function. We now demon- 
strate how the concept of resultants from elimination theory can be used to obtain the polynomial whose 
zero set is the required algebraic function. 

Definition 4.7 (Resultant). Given a polynomial 

a(x) = ao + a\ x + . . . + a„_i a;"" 1 + a n x n 
of degree n with roots at, for i = 1, . . . , n and a polynomial 

b(x) = 6 + h x + . . . + 6 m _! a;" 1 " 1 + b m x m 
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Operation: Ll v ,L^ v i — > L uv 


Matlab Code 


3 = < 

uv 


J v ffl u L 2 UV = det(ul - C^), where 
2 if Lj lv = L uv , 
^ (Cft ® + (lei ® C^) otherwise. 


function Luv3 = LlplusL2(Luvl,Luv2,u) 
Cul = Luv2Cu(Luvl,u); 
if (Luvl == Luv2) 

Cu3 = 2*Cul; 
else 

Cu2 = Luv2Cu(Luv2,u); 
Cu3 = kron(Cul,eye(length(Cu2))) + .. 
+kron(eye(length(Cul)) ,Cu2) ; 

end 

Luv3 = det(u*eye(length(Cu3))-Cu3); 


c^= < 


i v B u Z£ v = det(til- C£), where 
k C^ = CS;®CS» otherwise. 


function Luv3 = LltimesL2(Luvl,Luv2,u) 
Cul = Luv2Cu(Luvl,u); 
if (Luvl == Luv2) 

Cu3 = Cu2; 
else 

Cu2 = Luv2Cu(Luv2,u); 
Cu3 = kron(Cul,Cu2); 
end 

Luv3 = det(u*eye(length(Cu3))-Cu3); 



Table 5: Formal and computational description of the EB U and Kl u operators acting on the bivariate poly- 
nomials L uv (m, v) and L uv (w, v) where CJJ^ and CJJ^ arc their corresponding companion matrices 
constructed as in Table [3] and <g> is the matrix Kronccker product. 



of degree m with roots f3j, for j = 1, . . . , m, the resultant is defined as 

n m 

Res x (a(x) , b(x)) = a™ b n m [] [J(^ - a,). 

i=l 3=1 

From a computational standpoint, the resultant can be directly computed from the coefficients of the 
polynomials itself. The computation involves the formation of the Sylvester matrix and exploiting an identity 
that relates the determinant of the Sylvester matrix to the resultant. 

Definition 4.8 (Sylvester matrix). Given polynomials a(x) and b(x) with degree n and m respectively and 
coefficients as in Definition ^. 7[ the Sylvester matrix is the [n + m) x (n + rn) matrix 



S(a,b) 



a n 

0"n—l O n 



b m 

6 r „_i b m 



a ai 
a 







b h 

b 



Proposition 4.9. The resultant of two polynomials a{x) and b(x) is related to the determinant of the 
Sylvester matrix by 

det(S(a, b)) = Res x (a(x) , b(x)) 

Proof. This identity can be proved using standard linear algebra arguments. A proof may be found in [2]. 
□ For our purpose, the utility of this definition is that the EB U and M u operations can be expressed in terms 
of resultants. Suppose we are given two bivariate polynomials L uv and £ uv . By using the definition of the 
resultant and treating the bivariate polynomials as polynomials in u whose coefficients are polynomials in v, 
we obtain the identities 



Ll v (t,v) 



Res u (L uv (t-u,v),Ll v (u,v)) 



(4.2) 
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T 




pv 


= u 2 v + u (1 — v) + v 2 


1 

u 2 


1 v v 1 

• 1 " 
1 -1 • 
1 ■ 




" -v 
1 -l + v 

V 




-u 

1 — It + u 




L 2 UV = u 2 (v 2 - 3v + 1) + u (1 + v) + v 2 


1 

u 


1 V V 

■ 1 " 
1 1 ■ 

. 1 -3 1 . 














ir — o + 1 
L ij- — 3^ + 1- 


u -f- 1 

3 ?i — U 
It 2 + 1 













1 


u 2 


^ v 4 




o v V 8 














1 




■ 2 


-6 11 - 


10 18 


-8 1 " 






7-1 

-^uv 








u 

2 

II 

u 3 

u 4 


2 
5 
4 
1 


■ 2 

■ 1 


-8 4 
-4 2 
















1 


v v z 


v J 


v 4 


u° «° v' 


8 y 

V V 










1 














1 -6 


11 -6 1 _ 




7-1 

-^uv 




li 

u 3 
u 4 


-1 
1 


■ 1 

■ 1 


-1 


10 - 


1 3 

6 7-2 


-3 1 











7-1 


ffl v A™ 




r 2 




AU 


1 


1 


■ • • 1 " 




1 
u 


1 




?; 2 u 3 u 4 
■ 1 " 








.4 . . 




2 

U 






-2 1 




u 2 




■ 1 -4 ■ 




u 3 






. _4 . 




u 3 




-8 6 




u 4 


1 


i 


-9 3 ■ 




u* 


1 


-2 3 • • 




u 5 


2 - 


3 


7 ■ ■ 




u 5 


8 - 


-12 ■ ■ • 




u 6 


3 








u 6 


3 


2 ■ ■ ■ 




u 7 


4 




-1 • • 




u 7 


2 






u 8 


3 - 


1 


1 ■ ■ 




u 8 


-1 






u 9 
u 10 


2 
1 


3 







Tabic 6: Examples of EB and S3 operations on a pair of bivariate polynomials, and L\ 
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and 



L 3 uv (t 7 v)=Ll ® U L 2 U 



UV 



Res u {u D ^L 1 uv (t/u,v),Ll v (u,v)^ , 



(4.3) 



where D\ is the degree of L„ v with respect to u. By Proposition 14. 9[ evaluating the ES U and M u operations 
via the resultant formulation involves computing the determinant of the (D\ + D^) x (D^ + Sylvester 
matrix. When L„ v ^ L^ v , this results in a steep computational saving relative to the companion matrix 
based formulation in Table [5] which involves computing the determinant of a x (£)JflJ) matrix. Fast 

algorithms for computing the resultant exploit this and other properties of the Sylvester matrix formulation. 
In Maple , the computation L^ v = L^ v EB U L^ v may be performed using the command: 

Luv3 = subs (t=u, resultant (subs (u=t-u,Luvl) ,Luv2 ,u) ) ; 

The computation Lj^ v = L* v M u L^ v can be performed via the sequence of commands: 
Dul = degree (Luvl ,u) ; 

Luv3 = subs (t=u .result ant (simplify (u~Dul*subs (u=t/u,Luvl) ) ,Luv2 ,u) ) ; 

When L* v = L„ v , however, the EB U and £3 U operations are best performed using the companion matrix 
formulation in Tabled The software implementation of the operations in Table [5] in [22] uses the companion 
matrix formulation when L\ v = L^ v and the resultant formulation otherwise. 

Thus far we have established our ability to encode algebraic distribution as solutions of bivariate polynomial 
equations and to manipulate the solutions. This sets the stage for defining the class of "algebraic" random 
matrices next. 



We are interested in identifying canonical random matrix operations for which the limiting eigenvalue 
distribution of the resulting matrix is an algebraic distribution. This is equivalent to identifying operations 
for which the transformations in the random matrices can be mapped into transformations of the bivariate 
polynomial that encodes the limiting eigenvalue distribution function. This motivates the construction of 
the class of "algebraic" random matrices which we shall define next. 

The practical utility of this definition, which will become apparent in Section l6l and [PCTl can be succinctly 
summarized: if a random matrix is shown to be algebraic then its limiting eigenvalue density function can 
be computed using a simple root-finding algorithm. Furthermore, if the moments exist, they will satisfy a 
finite depth linear recursion (see Theorem 18. 6[) with polynomial coefficients so that we will often be able to 
enumerate them efficiently in closed form. Algebraicity of a random matrix thus acts as a certificate of the 
computability of its limiting eigenvalue density function and the associated moments. In this chapter our 
objective is to specify the class of algebraic random matrices by its generators. 

5.1 Preliminaries 

Let Ajv, for N = 1, 2, . . . be a sequence of N x N random matrices with real eigenvalues. Let F An 
denote the e.d.f., as in (|1.5|) . Suppose F An (x) converges almost surely (or in probability), for every x, to 
F A (x) as N — ► oo, then we say that Ajy *— > A. We denote the associated (non-random) limiting probability 
density function by Ja{x). 

Notation 5.1 (Mode of convergence of the empirical distribution function). When necessary we highlight 
the mode of convergence of the underlying distribution function thus: if An i— ^+ A then it is shorthand 
for the statement that the empirical distribution function of An converges almost surely to the distribution 
function F A ; likewise Ajy i— ^— * A is shorthand for the statement that the empirical distribution function of 
An converges in probability to the distribution function F A . When the distinction is not made then almost 
sure convergence is assumed. 



5. Class of algebraic random matrices 
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Remark 5.2. The element A above is not to be interpreted as a matrix. There is no convergence in the 
sense of an oo x oo matrix. The notation An <— — * A is shorthand for describing the convergence of the 
associated distribution functions and not of the matrix itself. We think of A as being an (abstract) element 
of a probability space with distribution function F A and associated density function f a ■ 

Definition 5.3 (Atomic random matrix). If f a G Vatom then we say that An is an atomic random matrix. 
We represent this as An i— * A £ M a tom where M a tom denotes the class of atomic random matrices. 

Definition 5.4 (Algebraic random matrix). If f a £ Taig then we say that An is an algebraically char- 
acterizable random matrix (often suppressing the word characterizable for brevity). We represent this as 
An 1 — ► A € Aiaig where M a i g denotes the class of algebraic random matrices. Note that, by definition, 

Matom dMalg- 

5.2 Key idea used in proving algebraicity preserving nature of a random matrix 
transformation 

The ability to describe the class of algebraic random matrices and the technique needed to compute the 
associated bivariatc polynomial is at the crux our investigation. In the theorems that follow, we accomplish 
the former by cataloguing random matrix operations that preserve algebraicity of the limiting distribution. 

Our proofs shall rely on exploiting the fact that some random matrix transformations, say An 1 — ► Bjv, 
can be most naturally expressed as transformations of L A Z i — ► L^ z ; others as Lf g i — ► Lf while some as 
L A y i — > L B y . Hence, we manipulate the bivariate polynomials (using the transformations depicted in Figure 
[3]) to the form needed to apply the appropriate operational law, which we derive as part of the proof, and then 
reverse the transformations to obtain the bivariate polynomial L^ z . Once we have derived the operational 
law for computing L^ z from L A Z , we have established the algebraicity of the limiting eigenvalue distribution 
of Bjv and we are done. Readers interested in the operational law may skip directly to Section [6?] 

The following property of the convergence of distributions will be invaluable in the proofs that follow . 

Proposition 5.5 (Continuous mapping theorem). Let An 1 — ► A. Let fA and S A denote the corresponding 
limiting density function and the atomic component of the support, respectively. Consider the mapping 
y = h(x) continuous everywhere on the real line except on the set of its discontinuities denoted by T>h- If 
Dh n S A = then B^ = /i(Ajv) 1 — ► B. The associated non-random distribution function, F B is given by 
F B (y) — F A (/i^ _1 ^(y)). The associated probability density function is its distributional derivative. 

Proof. This is a restatement of continuous mapping theorem which follows from well-known facts about 
the convergence of distributions [7] . □ 

5.3 Deterministic operations 

We first consider some simple deterministic transformations on an algebraic random matrix An that 
produce an algebraic random matrix B^. 



Theorem 5.6. Let An ^»4£ Maig and p, q, r, and s be real-valued scalars. Then, 

B N = (pA N + ql N )/(rA N + sI N ) i-> B e M alg , 
provided fA does not contain an atom at —s/r and r, s are not zero simultaneously. 



Proof. Here we have h(x) = (p x + r) /(q x + s) which is continuous everywhere except at x = —s/r for 
s and r not simultaneously zero. From Proposition 15. 5i unless /a(^) has an atomic component at —s/r, 
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B 



B. The Sticltjcs transform of F B can be expressed as 



m B (z) = E y 
Equation (|5.1[) can be rewritten as 

m B (z) = 



1 






= E X 


y - z_ 





r x + s 



p x + q — z(r x + s) 



rx + s 



(p — rz)x + (q — sz) 
With some algebraic manipulations, we can rewrite (|5.2p as 



dF A {x) = f ^±± F dF A {x 



(5.1) 



(5.2) 



m B (z) = f3 z 



r^±ldF A (x)=P z (r 



-dF A (x) + s 



1 



x + a 



-dF A (x) 



(3 z \r dF A {x)-ra, 



1 



(5.3) 



x + a 



-dF A (x) + s 



-dF A {x) 



where /3 Z = l/(p— r z) and a z = (q — sz)/(p— r z). Using the definition of the Stieltjes transform and the 
identity J dF A (x) — 1, we can express (tib(z) in (|5.3H in terms of m^z) as 



m B {z) = P z r + ((3 z s- (3ra z ) m A {~a z ). 
Equation (|5.4|) can, equivalently, be rewritten as 

m A (~a ) - m a( z ) ~ Pz r 

[3 z s - z r a z 

Equation (|5.5|) can be expressed as an operational law on L A Z as 



L^ z (m, z) = L A z ((m - /3 z r)/(/3 z s - f3 z ra z ), -a z ). 



(5.4) 



(5.5) 



(5.6) 



Since L A Z exists, we can obtain L^ lz by applying the transformation in (|5.6[) . and clearing the denominator 
to obtain the irreducible bivariate polynomial consistent with Remark 13. 31 Since L^ z exists, this proves that 
!b € Paig and Bjy i-> B e M a \ g . □ 

Appropriate substitutions for the scalars p, q, r and s in Theorem 1 5 . 61 leads to the following Corollary. 



Corollary 5.7. Let Am h4c M. a ig and let a be a real-valued scalar. Then, 

1. Bjv = A^y 1 \—>B^z Maig, provided f a does not contain at atom at 0, 

2. B N = a A N i-> B G M ; S) 

5. Ba, = A w ialjv i-» B e 7W ais . 



Theorem 5.8. Let X„jv be an n x N matrix. If Ajv = X„jvX n ^nic i/ien 

Bjv = X^ ^rXn^ I > -B £ A'laZg . 



Proof. Here X„jv is an ?i x JV matrix, so that A„ and Bjy are nxn and N x N sized matrices respectively. 
Let = n/A. When cat < 1, B^v will have N — n eigenvalues of magnitude zero while the remaining n 
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eigenvalues will be identically equal to the eigenvalues of A n . Thus, the e.d.f. of Bjy is related to the c.d.f. 
of A„ as 



N — n n „a , , 

= (1 - c N ) I [0)Oo) + cat F A " (x). 



(5.7) 



where I[o,oo) i s the indicator function that is equal to 1 when x > and is equal to zero otherwise. 

Similarly, when cn > 1, A„ will have n — N eigenvalues of magnitude zero while the remaining N 
eigenvalues will be identically equal to the eigenvalues of B^. Thus the e.d.f. of A„ is related to the e.d.f. 
of B 7v as 



F n (x) = I [0iOO ) + -F B »(x) 

n n 

= (i-—)i [Qoc) + — F Bn (x). 

\ C N C N 



(5.8) 



Equation (|5.8|) is (|5.7|l rearranged; so we do not need to differentiate between the case when cn < 1 and 

CN > 1. 

Thus, as n, N — > oo with cat = n/N — > c, if F A " converges to a non-random d.f. F A , then f Bjv will also 
converge to a non-random d.f. F B related to F A by 



F B (x) = (l-c)I [0iOo) +cF A (x). 



(5.9) 



From (|5.9p . it is evident that the Stieltjes transform of the limiting distribution functions F A and F B are 
related as 



m A (z) = - [ 1 - - ) - + -m B (z). 
\ c j z c 



(5.10) 



Rearranging the terms on either side of (|5.10j) allows us to express 7tib(z) in terms of m,A{z) as 

1 - c 



m B (z) 



cm A (z). 



(5.11) 



Equation (|5. 1 1[) can be expressed as an operational law on L A Z as 



L^ z (m, z) — L A Z ( - 


-H 




i \ 

— m, z . 








c J 



Given L A , we can obtain L^ z by using (|5.12D . Hence B^r | — *• B £ M a i g . □ 



(5.12) 



Theorem 5.9. Let A w i-> A G A1 a ; 9 . T/ien 



Bat = (Ajv) HBe^j, 



Proof. Here we have ft.(a;) = x 2 which is continuous everywhere. From Proposition 15.51 B^ 
Stieltjes transform of F B can be expressed as 



B. The 



m B (z) = E Y 



E 



x 



(5.13) 
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Equation (|5. 13|) can be rewritten as 



mB(z) 



2y / z J x — \fz 



dF A (x) - 



2y/~Z J X + sJ~Z 



dF (x) 



Equation (|5.14[) leads to the operational law 



im Z (™, z) = L A z (2my/z, y/z) ffl m L A z (-2m^/z, y/z). 



(5.14) 
(5.15) 

(5.16) 



Given L A Z , we can obtain L^ z by using (|5.16|) . This proves that Bjv^Be M a \ & . □ 



Theorem 5.10. Let A„ ^Ae M a i g and Bjv i— > B e M a i g - Then, 

C M = diag(A„, B N ) ^Ce A4 a ; s , 
where M = n + N and n/N — > c > as n, N oo. 



Proof. Let Cat be an N x N block diagonal matrix formed from the n x n matrix A n and the M x M 
matrix Bjv/- Let cn = n/N . The c.d.f. of Cjy is given by 

pC N 



cn 



F A " + (1 - c n )F Bm . 



Let n, N — > oo and cat = n/A^ — > c. If F An and i* 13 " converge in distribution almost surely (or in probability) 
to non-random d.f.'s F A and F B respectively, then F Cn will also converge in distribution almost surely (or 
in probability) to a non-random distribution function F c given by 



F c (x) = cF A (x) + (1 - c) F B {x). 



(5.17) 



The Stieltjes transform of the distribution function F c can hence be written in terms of the Stieltjes trans- 
forms of the distribution functions F A and F B as 

mc(z) = cm,A{z) + (1 - c) m B {z) (5.18) 

Equation (|5.18p can be expressed as an operational law on the bivariate polynomial L A z (m, z) as 



L mr , I — , z 



L a I — z 



(5.19) 



Given L A Z and and the definition of the EB m operator in Section L^ lz is a polynomial which can be 
constructed explicitly. This proves that C» ^ C e M a i g . □ 



Theorem 5.11. If A n = diag(Bjv, a I n _jv) and a is a real valued scalar. Then, 

BjV^5e Maig, 

as n, N — ► oo with cn — n/N — > c, 



Proof. Assume that as n, TV — » oo, cat = n/A~ — > c. As we did in the proof of Theorem 15. 101 we can show 
that the Stieltjes transform m,A(z) can be expressed in terms of msiz) as 



m A (z) 



a — z 



■ m B (z). 



(5.20) 
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This allows us to express L^ z (m, z) in terms of L^ z (m, z) using the relationship in (|5.20[) as 



Lg z (m,z)=Li z . 



1 1 

1 TO, Z 

a — z c 



We can hence obtain L^ z from using (|5.2ip . This proves that Bjv i— > B £ *M a ig- D 



(5.21) 



Corollary 5.12. Lei Ajv nie .Maig. ?7ien 

B w = diag(A„,aIj V -, l ) «B£ 7W a ; s , 

for n/N — > c > as n, N ^ oo. 



PROOF. This follows directly from Theorem 15. 101 □ 
5.4 Gaussian-like operations 

We now consider some simple stochastic transformations that "blur" the eigenvalues of A^r by injecting 
additional randomness. We show that canonical operations involving an algebraic random matrix Ajy and 
Gaussian-like and Wishart-like random matrices (defined next) produce an algebraic random matrix Bjv- 



Definition 5.13 (Gaussian-like random matrix). Let ~Yn,l be an N x L matrix with independent, identically 
distributed (i.i.d.) elements having zero mean, unit variance and bounded higher order moments. We label 
the matrix Gn.l = n.l as a Gaussian-like random matrix. 

We can sample a Gaussian-like random matrix in Matlab as 
G = sign(randn(N,L))/sqrt(L) ; 

Gaussian-like matrices are labelled thus because they exhibit the same limiting behavior in the N — > oo limit 
as "pure" Gaussian matrices which may be sampled in Matlab as 

G = randn(N,L)/sqrt(L) ; 

Definition 5.14 (Wishart-like random matrix). Let Gn,l be a Gaussian-like random matrix. We label the 
matrix Wn = Gn,l X G' N L as a Wishart-like random matrix. Let cn = N/L. We denote a Wishart-like 
random matrix thus formed by Wjv(cjv'). 



Remark 5.15 (Algebraicity of Wishart-like random matrices). The limiting eigenvalue distribution of the 
Wishart-like random matrix has the Marcenko-Pastur density which is an algebraic density since LY[ exists 



(see Table 1(b)). 



Assume that Gn,l is an N x L Gaussian-like random matrix. Let A^^—^A be an 

a.s. 



Proposition 5.16 

N x N symmetric /Hermitian random matrix and T^hAT be an L x L diagonal atomic random matrix 
respectively. //Gjv,l, An andTL are independent then B jy = Ajy + G N ^TlGjv.l^— — > B , ascL = N/L—*c 
forN,L — * oo,. The Stieltjes transform mg(z) of the unique distribution function F B is satisfies the equation 



mjfz) = tua z — c 



xdF T (x) 
1 + xmsiz) 



(5.22) 
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Proof. This result may be found in Marcenko-Pastur [18] and Silverstein [26]. □ 

We can reformulate Proposition 15 . 1 61 to obtain the following result on algebraic random matrices. 



Theorem 5.17. Let An, Gn,l and Tl be defined as in Provosition \5.16[ Then 

Bat = A w + G LN T L G L ^ N B e M a i g , 
as cl = N / L — > c for N, L — ► oo. 



Proof. Let T^ be an atomic matrix with d atomic masses of weight pi and magnitude A^ for i = 1, 2, . . . , d. 
From Proposition 1 5 . 1 61 ms(z) can be written in terms of m^fz) as 

m B (z) =m A \ z-cV . , ( \ )- ( 5 - 23 ) 

where we have substituted F T (x) = Y?,i=iPi^[K,°°) m to (|5.22p with ^2, i Pi = 1. 

Equation (|5 . 23[) can be expressed as an operational law on the bivariate polynomial L^ z as 



£ mz( m . Z ) = 4nz( m > z ~ «m)- 



(5.24) 



where a m = cj^ j=1 Aj/(1 + Ai to). This proves that B^v B € .Maig- □ 



Proposition 5.18. Assume that Wn(cn) is an N x N Wishart-like random matrix. Let An A be 
an N x N random Hermitian non-negative definite matrix. If Wjv(cjv) and An are independent, then 
Bjy = Ajy x Wjv(cAr) i? as cat — > c. The Stieltjes transform m B (z) 0/ the unique distribution function 
F B satisfies 



m B (z) = I dF {X) . (5.25) 

J \1 — c — c z m B (z)\x — z 

Proof. This result may be found in Bai and Silverstein [4,26]. □ 

We can reformulate Proposition 15 . 1 81 to obtain the following result on algebraic random matrices. 



Theorem 5.19. Let An and Wjv(c/v) satisfy the hypothesis of Provosition [5.18[ Then, 

B N = A N X W N {C N ) B G Malg, 

as c/v — * c. 



Proof. By rearranging the terms in the numerator and denominator, (|5 . 25[) can be rewritten as 

1 f dF A (x) 
m B (z) = j-r / V (5.26 

l-c-czm B (z) J x- x _ c _ czmB[z) 

Let a m ,z = 1 — c — czm B (z) so that (|5.26p can be rewritten as 

1 f dF A (x) 

m B {z) = / —777 7- (5-27) 
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Wc can express wb(z) in (|5.27p in terms of rriA(z) as 

1 



m B (z) 



m A (z/a„ hz ). 



Equation (|5.28j) can be rewritten as 

mA{z/a m , z ) = a m , 2 m B (z). 
Equation (|5.29j) can be expressed as an operational law on the bivariate polynomial L^ z as 



z) = i mz (a„ liZ to, z/a m<z ). 



(5.28) 
(5.29) 
(5.30) 



This proves that Bjy >— >-B £ M a \ g . □ 

Proposition 5.20. Assume that Gn.l is an N x L Gaussian-like random matrix. Let An^—^A be an 

1/2 

N x N symmetric/ Hermitian random matrix independent of Gn,l, A at. Let A^ denote an N x L matrix. 

X / 2 X / 2 / 3, s 

If s is a positive real-valued scalar then B N = (A N + y^s Gn,l )(A^ + ^G N , L ) ^B, as cl = N/L — > c 
/or N, L — > oo. T7ie Stieltjes transform, tob(z) of i/ie unique distribution function F B satisfies the equation 

dF A {x) 



m B {z) 



:{1 + scra B (z)} 



l+s c 



(T) + fl (C - 1) 



(5.31) 



Proof. This result is found in Dozier and Silverstein [12]. □ 

Wc can reformulate Proposition 15 . 201 to obtain the following result on algebraic random matrices. 



Theorem 5.21. Assume An, Gn,l and s satisfy the hypothesis of Provosition [5.2(A Then 

B N = (A# 2 + yfs G A r, i )(A^ /2 + yfs G N , L )' Sfie M alg , 
as Cl = N/L — > c for N, L — > oo . 



Proof. By rearranging the terms in the numerator and denominator, (|5.31[) can be rewritten as 

(z)= f {l + scm B {z)}dF A {x) 

Let a m = l + s cm B {z) and f3 m = {l + s cm B (z)}(z {l + s cm B (z)} + (c— 1) s), so that fj = af n z + a rn s(c— 1). 
Equation (|5.32|) can hence be rewritten as 



(5.32) 



m B (z) = QVi 



dF A (:r) 



(5.33) 



Using the definition of the Stieltjes transform in (|2.ip . we can express m B (z) in (|5.33|) in terms of toa(^) as 

tob(z) = a m m A (f3 m ) 

= a m TOa(o4 2 + a m (c - 1)8). 
Equation (|5.34|) can, equivalently, be rewritten as 

1 



(5.34) 



m A (a m z + a m {c - l)s) 



-m B {z). 



Equation (|5.35p can be expressed as an operational law on the bivariate polynomial L mz as 



L mz (m, z) = L mz (m/a m , a 2 z + a m s(c - 1)). 



(5.35) 



(5.36) 



This proves that B N B e M 



alg- 
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5.5 Sums and products 

Proposition 5.22. Let An 1 — ^ A and Hpft~-^B be N x N symmetric /Hermitian random matrices. Let 
Qjv be a Haar distributed unitary /orthogonal matrix independent of An and Bjv. Then Cn = An + 
QatBtvQat h ^ C The associated distribution function F c is the unique distribution function whose R 
transform satisfies 

rc(g) =r A (g)+r B (g)- (5.37) 
Proof. This result was obtained by Voiculescu in [34]. □ 

We can reformulate Proposition 15 . 221 to obtain the following result on algebraic random matrices. 



Theorem 5.23. Assume that An, Bat and Qn satisfy the hypothesis of Proposition [5. 22i Then, 

C N = A N + QjvBjvQat C G M a lg 



PROOF. Equation (|5.37[) can be expressed as an operational law on the bivariate polynomials Lf g and Lf & as 

(5.38) 



L Y g EB r L r 



If I/ mz exists then so does L rg and vice- versa. This proves that Cn h ^ C G -M a ig- □ 

Proposition 5.24. Let An^-^A and B^r h ^ B be N x N symmetric /Hermitian random matrices. Let 
Qjv be a Haar distributed unitary/orthogonal matrix independent of An and Bat. Then Cn = An X 
QatBtvQjv h ^ C where Cn is defined only if Cn has real eigenvalues for every sequence An and B^r. The 
associated distribution function F c is the unique distribution function whose S transform satisfies 

sc(y) = sa(u)sb(v)- (5.39) 
Proof. This result was obtained by Voiculescu in [35,36]. □ 

We can reformulate Proposition 15 . 241 to obtain the following result on algebraic random matrices. 



Theorem 5.25. Assume that An, and Bat satisfy the hypothesis of Proposition \5.24\ Then 

C N = A N x QatBtvQat C € M a ig ■ 



Proof. Equation (|5.39[) can be expressed as an operational law on the bivariate polynomials Lf y and L* 

(5.40) 



as 



If I/ mz exists then so docs L sy and vice versa. This proves that Bat >— > fi £ M a i g . □ 

Definition 5.26 (Orthogonally/Unitarily invariant random matrix). If the joint distribution of the elements 
of a random matrix An is invariant under orthogonal/unitary transformations, it is referred to as an or- 
thogonally /unitarily invariant random matrix. 
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If Ajy (or Bj\r) or both are an orthogonally/unitarily invariant sequences of random matrices then Theorems 
15.231 and 15.251 can be stated more simply. 



Corollary 5.27. Let An i-^— > A G M a ig and~f$N — * B>——* M a ig be a orthogonally/unitarily invariant random 
matrix independent of An- Then, 

1. C N = A N + B N ^CeMaig 

2. C N = A N x Bat h^-> C G M a i g 

Here multiplication is defined only if Cn has real eigenvalues for every sequence An and Bat. 



When both the limiting eigenvalue distributions of An and Bjv have compact support, it is possible to 
strengthen the mode of convergence in Theorems 15.231 and 15.251 to almost surely [15]. We suspect that al- 
most sure convergence must hold when the distributions are not compactly supported; this remains an open 
problem. 

6. Operational laws on bivariate polynomials 

The key idea behind the definition of algebraic random matrices in Section [5T1 was that when the limiting 
eigenvalue distribution of a random matrix can be encoded by a bivariate polynomial, then for the broad class 
of random matrix operations identified in Section [5j algebraicity of the eigenvalue distribution is preserved 
under the transformation. 

These operational laws, the associated random matrix transformation and the symbolic Matlab code 
for the operational law are summarized in Tables 091 The remainder of this chapter discusses techniques for 
extracting the density function from the polynomial and the special structure in the moments that allows 
them to be efficiently enumerated using symbolic methods. 



o 



p A+gl 
r A+s I 



A" 1 

A + a I 
a A 



A 
al 



B 
al 



A + G TG 



Operation 



B (m,z) 



MATLAB code 



Deterministic Transformations 



"Mobius' 



"Invert" 

'Translate" 
"Scale" 



' 'Projection/ 
Transpose" 

'Augmentation ' 



A x W(c) 



x 



"Add 
Atomic Wishart ' 



"Multiply 
Wishart' ' 



"Grammian' 



r A ( rn-(3 z r 

where a z = (q — s z)/(p — r z), 
and (3 Z = l/(p -rz). 



L ™ {am, 



Size of A 
Size of B 



c> 1 



T A 



m 

z c 



Size of A 
Size of B 



c < 1 



function LmzB - mobiusA(LmzA,p,q,r,s) 
syms m z 

alpha - ((q-s*z)/(p-r*z);beta=l/(p-r*z); 

temp_pol = subs(LmzA,z, -alpha); 

temp_pol - subs(temp_pol,m,((m/beta)-r)/(s-r*alpha)); 

LmzB - irreducLuv(temp_pol,m,z); 



function LmzB - invA(LmzA) 
LmzB - mobiusA(LmzA,0,l>l>0); 

function LmzB - shift A(LmzA, alpha) 
LmzB = mobiusA(LmzA,l,alpha,0,l); 

function LmzB - scaleA(LmzA) 
LmzB = mobiusA(LmzA,alpha,0,0,O; 



function LmzB = projectA(LmzA,c, alpha) 
syms m z 

mb = (l-(l/c))*(l/(alpha-z))+m/c; 
temp_pol = subs(LmzA,m,mb); 
LmzB - irreducLuv(temp_pol,m,z); 

function LmzB - augmentA(LmzA,c, alpha) 
syms m z 

mb = (l-(l/c))*(l/(alpha-z))+m/c; 
temp_pol = subs(LmzA,m,mb); 
LmzB - irreducLuv(temp_pol,m,z); 



Stochastic Transformations 



where a m = cYf i=1 
with J2i Pi ~ 



where a miZ = (1 — c — c z m). 



L mz ( — i + S(C - 1) ) , 

where a m = 1 + s cm. 



function LmzB = AplusWish(LmzA,c,p,lambda) 
syms m z 

alpha = z-c*sum(p.*(lambdaV(l-i-lambda*m))); 
temp_pol - subs(LmzA,z,z-alpha); 
LmzB - irreducLuv(temp_pol,m,z); 



function LmzB - AtimesWish(LmzA,c) 
syms mzzl 

alpha - (l-c-c*zl*m); temp_pol - subs(LmzA,m,m*alpha); 
temp_pol = subs(temp_pol,z,zl /alpha); 
temp_pol - subs(temp_pol,zl,z); % Replace dummy variable 
LmzB - irreducLuv(temp_pol,m,z); 



function LmzB - AgramG(LmzA,c,s) 
syms m z 

alpha - (l+s*c*m); beta - alpha*(z*alplia+s*(c-l)); 
temp_pol - subs(subs(LmzA,m,m/alpha),z,beta); 
LmzB - irreducLuv(temp_pol,m,z); 



V. 
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(a) L£ z ,— Lg z for A _ B = A 2 . 



Operational Law 


Matlab Code 


T A 

/ \ 

\ / 

fflm 
1 


function LmzB = squareA (LmzA) 
syms m z 

Lmzl = subs (LmzA, z, sqrt (z) ) ; 
Lmzl = subs (Lmzl ,m, 2*m*sqrt (z) ) ; 
Lmz2 = subs (LmzA, z , -sqrt (z) ) ; 
Lmz2 = subs (Lmz2 ,m, -2*m*sqrt (z) ) ; 

LmzB = LlplusL2(Lmzl ,Lmz2,m) ; 
LmzB = irreducLuv(LmzB,m,z) ; 



(b) L&x, L^ z i — > Lg z for A, B i > C = diag(A, B) where Size of A/ Size of C -> c. 



Operational Law 


Matlab Code 


t A jB 

mz mz 

I I 

\ / 

fflm 
1 

^mz 


function LmzC = AblockB (LmzA, LmzB, c) 
syms m z mu 

LmzAl = subs(LmzA,m,m/c) ; 
LmzBl = subs (LmzB, m, m/ ( 1-c) ) ; 

LmzC = LlplusL2(LmzAl,LmzBl,m) ; 
LmzC = irreducLuv(LmzC,m,z) ; 



Table 8: Operational laws on the bivariate polynomial encodings for some deterministic random matrix 
transformations. The operations EB U and M n are defined in Table 
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(a) Z4 z , L* z h_> Lg z for A, B i > C = A + QBQ' . 



Operational Law 


Matlab Code 


mz mz 

1 1 

tA tB 

rg rg 

\ / 

ffl r 

4 

1 

^mz 


function LmzC = AplusB(LmzA,LmzB) 
syms m z r g 

LrgA = Lmz2Lrg(LmzA) ; 
LrgB = Lmz2Lrg(LmzB) ; 

LrgC = LlplusL2 (LrgA, LrgB, r) ; 

LmzC = Lrg2Lmz(LrgC) ; 


(b) Z^ z , Lg z .— . Lg z for A, B i — > C = A x QBQ' . 


Operational Law 


Matlab Code 


mz mz 

4 1 

tA tB 

^ s y ^ S y 

\ / 
1 

r c 

1 


function LmzC = AtimesB(LmzA,LmzB) 
syms m z s y 

LsyA = Lmz2Lsy (LmzA) ; 
LsyB = Lmz2Lsy (LmzB) ; 

LsyC = LltimesL2 (LsyA, LsyB, s) ; 

LmzC = Lsy2Lmz(LsyC) ; 



Table 9: Operational laws on the bivariate polynomial encodings for some canonical random matrix trans- 
formations. The operations EB U and Kl u are defined in Table [5] 
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7. Interpreting the solution curves of polynomial equations 

Consider a bivariate polynomial L mz . Let D m be the degree of L mz (m, z) with respect to m and lk{z), for 
k = 0, . . . , D m , be polynomials in z that are the coefficients of m k . For every z along the real axis, there are 
at most D m solutions to the polynomial equation L mz (m, z) = 0. The solutions of the bivariate polynomial 
equation L mz = define a locus of points (to, z) in C x C referred to as a complex algebraic curve. Since 
the limiting density is over R, we may focus on real values of z. 

For almost every z € R, there will be D m values of m. The exception consists of the singularities of 
L m (m, z). A singularity occurs at z — zq if: 

• There is a reduction in the degree of to at zq so that there are less than D m roots for z = zq. This 
occurs when ZD m (zo) = 0. Poles of L mz (m, z) occur if some of the m-solutions blow up to infinity. 

• There are multiple roots of L mz at zq so that some of the values of to coalesce. 

The singularities constitute the so-called exceptional set of L mz (m, z). Singularity analysis, in the context 
of algebraic functions, is a well studied problem [14] from which we know that the singularities of L^ z (m, z) 
are constrained to be branch points. 

A branch of the algebraic curve L mz (m, z) = is the choice of a locally analytic function rrij(z) defined 
outside the exceptional set of L^ z (m,z) together with a connected region of the C x R plane throughout 
which this particular choice rrij(z) is analytic. These properties of singularities and branches of algebraic 
curve are helpful in determining the atomic and non-atomic component of the encoded probability density 
from L mz . We note that, as yet, we do not have a fully automated algorithm for extracting the limiting 
density function from the bivariate polynomial. Development of efficient computational algorithms that 
exploit the algebraic properties of the solution curve would be of great benefit to the community. 

7.1 The atomic component 

If there are any atomic components in the limiting density function, they will necessarily manifest them- 
selves as poles of L mz (m, z). This follows from the definition of the Stieltjes transform in (|2.ip . As mentioned 
in the discussion on the singularities of algebraic curves, the poles arc located at the roots of l-£> m (z). These 
may be computed in Maple using the sequence of commands: 

> Dm := degree (LmzA , m) ; 

> lDmz := coef f ( Lmz A, m, Dm) ; 

> poles := solve (lDmz=0, z) ; 

We can then compute the Puiseux expansion about each of the poles at z — zq. This can be computed 
in Maple using the algcurves package as: 

> with(algcurves) : 

> puiseux (Lmz ,z=pole ,m, 1) ; 

For the pole at z = zq, we inspect the Puiseux expansions for branches with leading term 1/(zq — z). 
An atomic component in the limiting spectrum occurs if and only if the coefficient of such a branch is non- 
negative and not greater than one. This constraint ensures that the branch is associated with the Stieltjes 
transform of a valid probability distribution function. 

Of course, as is often the case with algebraic curves, pathological cases can be easily constructed. For 
example, more than one branch of the Puiseux expansion might correspond to a candidate atomic compo- 
nent, i.e., the coefficients arc non-negative and not greater than one. In our experimentation, whenever 
this has happened it has been possible to eliminate the spurious branch by matrix theoretic arguments. 
Demonstrating this rigorously using analytical arguments remains an open problem. 

Sometimes it is possible to encounter a double pole at z = z corresponding to two admissible weights. 
In such cases, empirical evidence suggests that the branch with the largest coefficient (less than one) is the 
"right" Puiseux expansion though we have no theoretical justification for this choice. 
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7.2 The non-atomic component 

The probability density function can be recovered from the Stieltjes transform by applying the inversion 
formula in (|2.4p . Since the Stieltjes transform is encoded in the bivariate polynomial L mz , we accomplish 
this by first computing all D m roots along z € R (except at poles or singularities). There will be D m roots 
of which one solution curve will be the "correct" solution , i.e., the non-atomic component of the desired 
density function is the imaginary part of the correct solution normalized by n. In Matlab , the D m roots 
can be computed using the sequence of commands: 

Lmz_roots = [] ; 

x_range = [x_start :x_step:x_end] ; 
for x = x_range 

Lmz_roots_unnorm = roots (sym2poly (subs(Lmz ,z ,x) ) ) ; 

Lmz_roots = [Lmz_roots; 

real (Lmz_roots_unnorm) + i*imag(Lmz_roots_unnorm)/pi] ; 

end 

The density of the limiting eigenvalue distribution function can be, gencrically, be expressed in closed 
form when D m = 2. When using root-finding algorithms, for D m = 2, 3, the correct solution can often be 
easily identified; the imaginary branch will always appear with its complex conjugate. The density is just 
the scaled (by 1/tt) positive imaginary component. 

When D m > 4, except when L mz is bi-quadratic for D m = 4, there is no choice but to manually isolate 
the correct solution among the numerically computed D m roots of the polynomial Lm Z m, z) at each z = zq. 
The class of algebraic random matrices whose eigenvalue density function can be expressed in closed form is 
thus a much smaller subset of the class of algebraic random matrices. When the underlying density function 
is compactly supported, the boundary points will be singularities of the algebraic curve. 

In particular, when the probability density function is compactly supported and the boundary points 
are not poles, they occur at points where some values of m coalesce. These points are the roots of the 
discriminant of L mz , computed in Maple as: 

> PossibleBoundaryPoints = solve(discrim(Lmz,m) ,z) ; 

We suspect that "nearly all" algebraic random matrices with compactly supported eigenvalue distribution 
will exhibit a square root type behavior near boundary points at which there are no poles. In the generic case, 
this will occur whenever the boundary points correspond to locations where two branches of the algebraic 
curve coalesce. 

For a class of random matrices that includes a subclass of algebraic random matrices, this has been 
established in [27]. This endpoint behavior has also been observed orthogonally /unitarily invariant random 
matrices whose distribution has the element-wise joint density function of the form 

/(A) = C N exp (-NTrV(A)) dA 

where V is an even degree polynomial with positive leading coefficient and dA is the Lebesgue measure 
on N x N symmetric/Hermitian matrices. In [9], it is shown that these random matrices have a limiting 
mean eigenvalue density in the N — > oo limit that is algebraic and compactly supported. The behavior at 
the endpoint typically vanishes like a square root, though higher order vanishing at endpoints is possible 
and a full classification is made in [10]. In [17] it is shown that square root vanishing is generic. A similar 
classification for the general class of algebraic random matrices remains an open problem. This problem is of 
interest because of the intimate connection between the endpoint behavior and the Tracy- Widom distribution. 
Specifically, we conjecture that "nearly all" algebraic random matrices with compactly supported eigenvalue 
distribution whose density function vanishes as the square root at the endpoints will, with appropriate 
re-centering and rescaling, exhibit Tracy- Widom fluctuations. 

Whether the encoded distribution is compactly supported or not, the —1/z behavior of the real part of 
Stieltjes transform (the principal value) as z — > ±oo helps isolate the correct solution. In our experience, 
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while multiple solution curves might exhibit this behavior, invariably only one solution will have an imagi- 
nary branch that, when normalized, will correspond to a valid probability density. Why this always appears 
to be the case for the operational laws described is a bit of a mystery to us. 



Example: Consider the Marcenko-Pastur density encoded by L mz given in Table 1(b) The Puiseux ex- 
pansion about the pole at z = (the only pole!), has coefficient (1 — 1/c) which corresponds to an atom 
only when c > 1 (as expected using a matrix theoretic argument). Finally, the branch points at (1 ± \/c) 2 
correspond to boundary points of the compactly supported probability density. Figure [4] plots the real and 
imaginary parts of the algebraic curve for c = 2. 



8. Enumerating the moments and free cumulants 

In principle, the moments generating function can be extracted from L MZ by a Puiseux expansion of the 
algebraic function p{z) about z = 0. When the moments of an algebraic probability distribution exist, there 
is additional structure in the moments and free cumulants that allows us to enumerate them efficiently. For 
an algebraic probability distribution, we conjecture that the moments of all order exist if and only if the 
distribution is compactly supported. 

Definition 8.1 (Rational generating function). Let R[[x]] denote the ring of formal power series (or gener- 
ating functions) in x with real coefficients. A formal power series (or generating function) v G R[[u]] is said 
to be rational if there exist polynomials in u, P{u) and Q{u), Q(0) ^ such that 

Q[u) 

Definition 8.2 (Algebraic generating function). Let denote the ring of formal power series (or gen- 

erating functions) in x with real coefficients. A formal power series (or generating function) v G K[[ti]] is 
said to be algebraic if there exist polynomials in u, Pq(u), . . . , Pd u (u), n °t a ^ identically zero, such that 

Pq(u) + Pi(u)v + ... + P Bv (u)v^ = 0. 

The degree of v is said to be D v . 

Definition 8.3 (D-finite generating function). Let v G R[[u]]. If there exist polynomials po(u), . . . ,pd(u), 
such that 

p d (u)vM +p d ^(u)v^ + . . .+ Pl (u)v^ + p (u) = 0, (8.1) 

where v^' = div/du^ . Then we say that v is a D-finite (short for differentiably finite) generating function 
(or power series). The generating function, v(u), is also referred to as a holonomic function. 

Definition 8.4 (P-recursive coefficients). Let a n for n > denote the coefficients of a D-finite series v. If 
there exist polynomials Pa, . . . , P e G K[n] with P e ^ 0, such that 

P e (n)a n+e + P e -i(n)a n+e -i + ■ ■ ■ + Po(n)a n = 0, 

for all n G N, then the coefficients a n are said to be P-recursive (short for polynomially recursive). 

Proposition 8.5. Let v G be an algebraic power series of degree D v . Then v is D-finite and satisfies 

an equation A8.1\) of order D v . 



PROOF. A proof appears in Stanley [30, pp.187]. □ 



The structure of the limiting moments and free cumulants associated with algebraic densities is described 
next. 
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^— ^— Extraneous solution 
- - - Correct solution 
















Y 
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z 

(a) Real component. The singularity at zero corresponds to an atom of weight 1/2. The branch points at 
(1 ± V2) 2 correspond to the boundary points of the region of support. 




Figure 4: The real and imaginary components of the algebraic curve defined by the equation L mz (m, z) = 0, 
where L mz = czm 2 — (1 — c — z) m + 1, which encodes the Marcenko-Pastur density. The curve is 
plotted for c = 2. The —1/z behavior of the real part of the "correct solution" as z — > 00 is the 
generic behavior exhibited by the real part of the Stieltjcs transform of a valid probability density 
function. 
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Theorem 8.6. If /a G Vaig, o,nd the moments exist, then the moment and free cumulant generating functions 
are algebraic power series. Moreover, both generating functions are D-finite and the coefficients are P- 
recursive. 



PROOF. If f a € "^aig, then exists. Hence L^ z and Lf g exist, so that Ha{z) and are algebraic power 

series. By Theorem 18.51 thev are D-finite; the moments and free cumulants are hence P-recursive. □ 



There are powerful symbolic tools available for enumerating the coefficients of algebraic power series. 
The Maple based package gfun is one such example [24]. From the bivariate polynomial L^ z , wc can obtain 
the scries expansion up to degree expansion_degree by using the commands: 

> with(gfun) : 

> MomentSeries = algeqtoseries(Lmyuz ,z ,myu, expansion_degree , 'pos_slopes ' ) ; 

The option pos_slopes computes only those branches tending to zero. Similarly, the free cumulants can 
be enumerated from L rg using the commands: 

> with(gfun) : 

> FreeCumulantSeries = algeqtoseries(Lrg,g,r,expansion_degree, 'pos_slopes' ) ; 

For computing expansions to a large order, it is best to work with the recurrence relation. For an algebraic 
power series v(u), the first number_of _terms coefficients can be computed from L uv using the sequence of 
commands: 

> with(gfun) : 

> deq := algeqtodif f eq(Luv,v(u) ) ; 

> rec := diff eqtorec(deq,v(u) ,a(n)) ; 

> p_generator := rectoproc(rec,a(n) ,list) : 

> p_generator (number_of _terms) ; 



Example: Consider the Marcenko-Pastur density encoded by the bivariate polynomials listed in Table 1(b) 
Using the above sequence of commands, we can enumerate the first five terms of its moment generating 
function as 

a{z) = 1 + z + (c + l)z 2 + (3c + c 2 + 1) z 3 + (6 c 2 + c 3 + 6 c + l) z 4 + O (z 5 ) . 



The moment generating function is a D-Finite power series and satisfies the second order differential equation 
-z + zc - 1 + (-Z - zc + 1) n (z) + (z 3 c 2 - 2 z 2 c - 2 z 3 c + z - 2 z 2 + z 3 ) -j-fx (z) = 0, 

with initial condition /n(0) = 1. The moments M n = a(n) themselves are P-recursive satisfying the finite 
depth recursion 

(-2 c + c + 1) na (n) + ((-2 - 2 c) n - 3 c - 3) a (n + 1) + (3 + n) a (n + 2) = 



with the initial conditions a (0) = 1 and a (I) = 1. The free cumulants can be analogously computed. 

What we find rather remarkable is that for algebraic random matrices, it is often possible to enumerate 
the moments in closed form even when the limiting density function cannot. The linear recurrence satisfied 
by the moments may be used to analyze their asymptotic growth. 

When using the sequence of commands described, sometimes more than one solution might emerge. In 
such cases, we have often found that one can identify the correct solution by checking for the positivity of 
even moments or the condition /i(0) = 1. More sophisticated arguments might be needed for pathological 
cases. It might involve verifying, using techniques such as those in [1], that the coefficients enumerated 
correspond to the moments a valid distribution function. 
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9. Computational free probability 

9.1 Moments of random matrices and asymptotic freeness 

Assume we know the eigenvalue distribution of two matrices A and B. In general, using that information 
alone, we cannot say much about the the eigenvalue distribution of the sum A + B of the matrices since 
eigenvalues of the sum of the matrices depend on the eigenvalues of A and the eigenvalues of B, and also 
on the relation between the eigenspaces of A and of B. However, if we pose this question in the context of 
N x TV-random matrices, then in many situations the answer becomes deterministic in the limit N — > oo. 
Free probability provides the analytical framework for characterizing this limiting behavior. 

Definition 9.1. Let A = (Ajv)jveN be a sequence of N x N -random matrices. We say that A has a limit 
eigenvalue distribution if the limit of all moments 

a n := lim E[tr(Afr)] (n G N) 

N^oo 

exists, where E denotes the expectation and tr the normalized trace. 

Using the language of limit eigenvalue distribution as in Definition 19.11 our question becomes: Given 
two random matrix ensembles of N x TV-random matrices, A = {K^Nen and B = (Bjv)jveN, with limit 
eigenvalue distribution, does their sum C = (Cn)ngn, with Cat = A-n + Bjv, have a limit eigenvalue 
distribution, and furthermore, can we calculate the limit moments a% of C out of the limit moments (a£)k>i 
of A and the limit moments (aj? )k>i of B in a deterministic way. It turns out that this is the case if the 
two ensembles are in generic position, and then the rule for calculating the limit moments of C are given by 
Voiculescu's concept of "freeness." 

Theorem 9.2 (Voiculescu [36]). Let A and B be two random matrix ensembles of N x N -random matrices, 
A = (Ajv)./vgN and B = (Bn)ngn, each of them with a limit eigenvalue distribution. Assume that A and 
B are independent (i.e., for each N G N, all entries of A^r are independent from all entries o/Bjvj, and 
that at least one of them is unitarily invariant (i.e., for each N, the joint distribution of the entries does 
not change if we conjugate the random matrix with an arbitrary unitary N x N matrix). Then A and B are 
asymptotically free in the sense of the following definition. 

Definition 9.3 (Voiculescu [33]). Two random matrix ensembles A = (A^at^n and B = (Btv)at£n with 
limit eigenvalue distributions are asymptotically free if we have for all p > 1 and all n(l), m(l), . . . , n(p), 
m(p) > 1 that 



lim E 



tr{(A 



t (i) 

N 



^(l) 1 ) 



(B™ (1) 



a£ (1) l)...(A"<*> 



"n(p) 1 ) 



m(p) 



1)} 







In essence, asymptotic freeness is actually a rule which allows to calculate all mixed moments in A and 
B, i.e., all expressions of the form 

lim £[tr(A n WB m WA n(2 >B m ( 2 ) ■ • ■ A"^B m ^)] 

N—>oc 

out of the limit moments of A and the limit moments of B. In particular, this means that all limit moments 
of A + B (which are sums of mixed moments) exist, thus A + B has a limit distribution, and are actually 
determined in terms of the limit moments of A and the limit moments of B. For more on free probability, 
including extensions to the setting where the moments do not exist, we refer the reader to [6,15,21,37]. 

We now clarify the connection between the operational law of a subclass of algebraic random matrices 
and the convolution operations of free probability. This will bring into sharp focus how the polynomial 
method constitutes a framework for computational free probability theory 
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Free additive convolution 


Ja+b 


= fA 


ffl/fl 


L rg 




EBr -^^ 


Free multiplicative convolution 


fAxB 


= Sa 




tAxB 


- ^sy 





Table 10: Implicit representation of the free convolution of two algebraic probability densities. 



Proposition 9.4. Let Ajyi — >A and Hn 1 -^B be two asymptotically free random matrix sequences as in 
Definition GO Then A N + B N ^A + B and Ajv x Bjv i — ► AB ( where the product is defined whenever 
An x Byv has real eigenvalues for every An and Bjv ) with the corresponding limit eigenvalue density 
functions, fA+B and f ab given by 

f A+B = f A ffl f B (9.1a) 

f AB = f A ^ f B (9.1b) 

where EB denotes free additive convolution and M denotes free multiplicative convolution. These convolution 
operations can be expressed in terms of the R and S transforms as described in Propositions \5.2B and \5^\ 
respectively. 

Proof. This result appears for density functions with compact support in [34,35]. It was later strengthened 
to the case of density functions with unbounded support. See [15] for additional details and references. □ 

In Theorems 15.231 and 15 . 251 we. in effect, showed that the free convolution of algebraic densities produces 
an algebraic density. This stated succinctly next. 



Corollary 9.5. Algebraic probability distributions form a semi-group under free additive convolution. 

Corollary 9.6. Algebraic distributions with positive semi-definite support form a semi-group under free 
multiplicative convolution. 



This establishes a framework for computational free probability theory by identifying the class of distri- 
butions for which the free convolution operations produce a "computable" distribution. 

9.2 Implicitly encoding the free convolution computations 

The computational framework established relics on being able to implicitly encode free convolution com- 
putations as a resultant computation on appropriate bivariate polynomials as in Table 1101 This leads to the 
obvious question: Are there other more effective ways to implicitly encode free convolution computations? 
The answer to this rhetorical question will bring into sharp focus the reason why the bivariate polynomial 
encoding at the heart of the polynomial method is indispensable for any symbolic computational implemen- 
tation of free convolution. First, we answer the analogous question about the most effective encoding for 
classical convolution computations. 

Recall that classical convolution can be expressed in terms of the Laplace transform of the distribu- 
tion function. In what follows, we assume that the distributions have finite moments^. Hence the Laplace 
transform can be written as a formal exponential moment generating function. Classical additive and mul- 
tiplicative convolution of two distributions produces a distribution whose exponential moment generating 
function equals the series (or Cauchy) product and the coefficient-wise (or Hadamard) product of the in- 
dividual exponential moment generating functions, respectively. Often, however, the Laplace transform of 



In the general case, tools from complex analysis can be used to extend the argument. 
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cither or both the individual distributions being convolved cannot be written in closed form. The next 
best thing to do then is to find an implicit way to encode the Laplace transform and to do the convolution 
computations via this representation. 

When this point of view is adopted, the task of identifying candidate encodings is reduced to finding the 
class of representations of the exponential generating function that remains closed under the Cauchy and 
Hadamard product. Clearly rational generating functions (see Definition 18. lj) satisfy this requirement. It 
is shown in Theorem 6.4.12 [30, pp.194], that D-finitc generating functions (see Definition I8.3() satisfy this 
requirement as well. 

Proposition 18.51 establishes that all algebraic generating functions (sec Definition 18. 2[) and by extension, 
rational generating functions, are also D-finite. However, not all D-finite generating functions are algebraic 
(see Exercise 6.1 [30, pp. 217] for a counter-example) so that algebraic generating functions do not satisfy 
the closure requirement. Furthermore, from Proposition 6.4.3 and Theorem 6.4.12 in [30], if the ordinary 
generating function is D-finite then so is the exponential generating function and vice versa. Thus D-finitc 
generating functions are the largest class of generating functions for which classical convolution computations 
can be performed via an implicit representation. 

In the context of developing a computational framework based on the chosen implicit representation, it 
is important to consider computability and algorithmic efficiency issues. The class of D-finitc functions is 
well-suited in that regard as well [24] so that we regard it as the most effective class of representations in 
which the classical convolution computations may be performed implicitly. 

However, this class is inadequate for performing free convolution computations implicitly. This is a 
consequence of the prominent role occupied in this theory by ordinary generating functions. Specifically, 
the ordinary formal R and S power series, are obtained from the ordinary moment generating function by 
functional inversion (or reversion) , and are the key ingredients of free additive and multiplicative convolution 
(see Propositions 19 .41 15 . 221 and 15.24)) . The task of identifying candidate encodings is thus reduced to finding 
the class of representations of the ordinary moment generating function that remains closed under addition, 
the Cauchy product, and reversion. D-finitc functions only satisfy the first two conditions and are hence 
unsuitable representations. 

Algebraic functions do, however, satisfy all three conditions. The algorithmic efficiency of computing the 
resultant (see Section 14. 2[) justifies our labelling of the bivariatc polynomial encoding as the most effective 
way of implicitly encoding free convolution computations. The candidacy of constructibly D-finite generating 
functions [5], which do not contain the class of D-finite functions but do contain the class of algebraic 
functions, merits further investigation since they are closed under reversion, addition and multiplication. 
Identifying classes of representations of generating functions for which both the classical and free convolution 
computations can be performed implicitly and effectively remains an important open problem. 

10. Applications 

We illustrate the use of the computational techniques developed in Section [6] with some examples. 
Documented MATLAB implementation of the polynomial method is available via the RMTool package [22] 
from http://www.mit.edu/~raj/rmtool/; the examples considered in this article, along with many more, 
appear there and in [20]. 

10.1 The Jacobi random matrix 

The Jacobi matrix ensemble is defined in terms of two independent Wishart matrices Wi(ci) and W2(c2) 
as J = (I + W 2 (c 2 ) W 1 _1 (ci)) _1 . The subscripts are not to be confused for the size of the matrices. Listing 
the computational steps needed to generate a realization of this ensemble, as in Table QTJ is the easiest way 
to identify the sequence of random matrix operations needed to obtain L^ z . 
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Transformation 


Numerical MATLAB code 


Symbolic MATLAB code 


Initialization 


% Pick n, cl, c2 
lNl=n/cl; lNz=n/cz; 


% Define symbolic variables 
syms m c z; 


A, - T 


A 1 fiirn | n Y~i l ■ 

rt. J. U VL I 11) 11 J j 


T iTn7 l — m*(l ?\ 1 - 

-LjIII/^-L 111 \± J 1 


A 2 =Wi(ci) x Aj 


Gl = randn(n,Nl)/sqrt(Nl); 
Wl = Gl*Gl'; 
A2 = Wl*Al; 


Lmz2 = AtimcsWish(Lmzl,cl); 


A 3 = Aj 1 


A3 = inv(A2); 


Lmz3 = invA(Lmz2); 


A 4 = W 2 (c 2 ) x A 3 


G2 = randn(n,N2)/sqrt(N2); 
W2 = G2*G2'; 
A4 = W2*A3; 


Lmz4 = AtimesWish(Lmz3,c2); 


A 5 = A 4 + I 


A5 = A4+I; 


Lmz5 = shiftA(Lmz4,l); 


A 6 =A- 1 


A6 = inv(A5); 


Lmz6 = invA(Lmz5); 



Table 11: Sequence of MATLAB commands for sampling the Jacobi ensemble. The functions used to gen- 
erate the corresponding bivariate polynomials symbolically are listed in Table [7] 



We first start off with Ai = I. The bivariate polynomial that encodes the Stieltjcs transform of its eigenvalue 
distribution function is given by 

Lmz(m,z) = (1 - z)m - 1. 



For A2 = Wi(ci) x Ai, we can use (|5.30[) to obtain the bivariate polynomial 

L 2 nz (m, z) = z Ci m 2 — (— ci — z + 1) m + 1. 

For A 3 = A^ 1 , fr om (|5.6|) . we obtain the bivariate polynomial 

L^ nz (m, 2) = z 2 c\Tn 2 + (ci z + z — 1) m + 1. 

For A4 = W"2(c2) x A3. We can use (|5.30[) to obtain the bivariate polynomial 

L^ z (m, z) — (ci z 2 + C2 z) m 2 + (ci z + 2 — 1 + C2) m + 1. 

For A5 = A4 + I, from (|5.6p . we obtain the bivariate polynomial 



L mz (m, 2) = ((2 - l) 2 ci + c-2 (2 — 1)) m + (ci (2 — 1) + 2 — 2 + c 2 ) m + 1. 



Finally, for J = A6 = A 5 , from (J5T6J) , we obtain the required bivariate polynomial 



Lmz(m, 2) = L mz (m, 2) = (c x 2 + 2 3 ci — 2ci 2 2 



3 . 2\ 2 

c 2 2 + c 2 2 ) m 



(10.1) 
(10.2) 

(10.3) 

(10.4) 

(10.5) 



+ (— 1 + 2 2 + ci — 3ci 2 + 2ci 2 2 + c 2 2 — 2 c 2 2 2 ) m — c 2 2 — ci + 2 + ci 2. (10.6) 

Using matrix theoretic arguments, it is clear that the random matrix ensembles A3, . . . Aq are defined only 
when ci < 1. There will be an atomic mass of weight (1 — I/C2) at 1 whenever C2 > 1. The non-atomic 
component of the distribution will have a region of support S n = (a_,a + ). The limiting density function 
for each of these ensembles can be expressed as 



fA t (x) 



y/{x- a-)(a+ - x) 
2tt/ 2 (x) 



for a_ < x < d-|_, 



(10.7) 



for i = 2,..., 6, where a_, a+ , where the polynomials hix) are listed in Table [121 The moments for 
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h{x) 






X C\ 


(i ± y/°rr 


A 3 


2 

X Cl 


i 

(1 -F v^T) 2 


A 4 


C±X 2 + C2X 


1 + ci + C2 - C1C2 ± 2-y/ci + C2 - C1C2 
(1-ci) 2 


A 5 


Cl {x- l) 2 +c 2 (x- 1) 


cf - cl + 2 + c 2 — cic 2 ± 2Vci + c 2 — cic 2 
(1-ci) 2 


A 6 


(ci x + x 3 ci — 2 ci x 2 — C2 x 3 + C2 i 2 ) 


(1-Cl) 2 

c 2 — cl + 2 + C2 — C1C2 =F 2-y/ci + C2 — C1C2 



Tabic 12: Parameters for determining the limiting eigenvalue density function using (|10.7p . 

1.6p 
1.4 - 
1.2 - 

1 - 

LL 08 
Q 

°- 0.6- 
0.4 
0.2 - 
0^ 

0.2 0.4 0.6 0.8 1 1.2 

x 

Figure 5: The limiting density (solid line), fA 6 {x), given by (| 10 . 7[) with c\ = 0.1 and c 2 = 0.625 is compared 
with the normalized histogram of the eigenvalues of a Jacobi matrix generated using the code in 
TableEDover 4000 Monte-Carlo trials with n = 100, N x = n/c x = 1000 and N 2 = n/c 2 = 160. 

the general case when c\ ^ c 2 can be enumerated using the techniques described; they will be quite messy. 
Instead, consider the special case when c\ — c 2 = c. Using the tools described, the first four terms of the 
moment series, n(z) = fij(z), can be computed directly from L^ z as 

r\ 1 ^f 1 -l 1 ^ ^f 3 j- 1 ^ 2 ^( X 2 ^ 3 1 3, 1\ a 
u(z) — — h - c H — \ z + I — cH — z + — c H c c H z 

2 \& 4J \W &) \32 16 128 16 J 

The moment generating function satisfies the differential equation 

-3 z + 2 + zc + (-6 z 2 + z 3 + 10 z + z 3 c 2 - 2 « 3 c - 4) fj. (z) 

+ (z 4 - 5 z 3 - 2 A + 8 z 2 + z 4 c 2 + 2 z 3 c - 4 2 - z 3 c 2 ) —ft («) = 0, 

with the initial condition /x(0) = 1. The moments a(n) = M n themselves are P-recursive and obtained by 




The polynomial method 



41 



the recursion 



(-2c + c 2 + 1 + (-2c + c 2 + 1) n)a{n) + ((-5 + 2 c - c 2 ) n - 11 4- 2 c - c 2 ) a (n + 1) 

4- (26 4- 8n) a (n 4- 2) + (-16 - 4n) a (n + 3) = 0, 

with the initial conditions a(0) = 1/2, a(l) = l/8c4-l/4, and a(2) = 3/16c4-l/8. We can similarly compute 
the recursion for the free cumulants, a(n) = K n+ \ 1 as 

nca (n) + (12 + 4 n) a (n + 2) = 0, 

with the initial conditions a(0) = 1/2, and a(l) = 1/8 c. 

10.2 Random compression of a matrix 



Theorem 10.1. Let Aat h- > A € "P a ig- £e< Qat be an N x N Haar unitary/orthogonal random matrix 
independent of An- Let B„ &e the upper n x n Wocfc o/QatAatQat- TTien 



as n/iV — > c for n, N 



CO. 



Proof. Let Pat be an AT x A 1 " projection matrix 



N 



Q 



N 



>N-n 



Q 



By definition, P^v is an atomic matrix so that Pjv — > P G A'laig as n/iV — > c for n, AT — > oo. Let Bjy = 
Ptv x Aat. By Corollarv l5.27l Bat — ► i? e jM a i g - Finally, from Theorem l5.1H we have that B„ ->B£ A4 a i g - 
□ 

The proof above provides a recipe for computing the bivariate polynomial explicitly as a function of 
L^ z and the compression factor c. For this particular application, however, one can use first principles [29] 
to directly obtain the relationship 

rB(g) = r A {cg), 



expressed in terms of the R transform. This translates into the operational law 



L ?A r ,g) = L t Jr,cg). 



(10.8) 



Example: Consider the atomic matrix An half of whose eigenvalues are equal to one while the remainder 
are equal to zero. Its eigenvalue distribution function is given by (|3.6p . From the bivariate polynomial, Lf E 
in Table 1(a) and (| 10 it can be show that the limiting eigenvalue distribution function of B„, constructed 



from Aat as in Theorem 1 10.H is encoded by the polynomial 



-2 cz + 2 cz) ra - (-2 c + 4 cz + 1 - 2 z) m - 2 c + 2, 



where c is the limiting compression factor. Poles occur at z = and z = 1. The leading terms of the Puiseux 
expansion of the two branches about the poles at z = zq are 



The polynomial method 



42 





Figure 6: The limiting eigenvalue density function (solid line) of the top 0.4iV x OAN block of a randomly 
rotated matrix is compared with the experimental histogram collected over 4000 trials with N = 
200. Half of the eigenvalues of the original matrix were equal to one while the remainder were 
equal to zero. 



It can be easily seen that when c > 1/2, the Puiscux expansion about the poles z = zo will correspond to an 
atom of weight wq = (2c — l)/2c. Thus the limiting eigenvalue distribution function has density 



tl\ / 2c-l n ^ 1 y/(x-a-)(a+-x) ( 2c-l \ 

fB{x) = max (-^.Oj 5(x) + - 2xc _ 2cx2 J [a _, a+1 +max ,0j 5{x 1), 



(10.9) 



where a± = 1/2 ± V— c 2 + c. Figure PT0.2I compares the theoretical prediction in (|10.9|) with a Monte-Carlo 
experiment for c = 0.4. From the associated bivariate polynomial 

L B = (-2 c + 2 cz) /i 2 + (z - 2 - 2 cz + 4 c) \i - 2 c + 2, 



we obtain two series expansions whose branches tend to zero. The first four terms of the series are given by 

l + ±z+±±Zz* + l±Zz* + 0(z% (10.10) 

and, 



C -l + £^l z _(^-l)(-2 + C ) z3 _( C -lH3 C -4) z , +o(A 



( 1(U1 ) 

respectively. Since c < 1, the series expansion in (|10.11|) can be eliminated since /x(0) := / dF B {x) = 1. 
Thus the coefficients of the series in (|10.10|) are the correct moments of the limiting eigenvalue distribution. 
A recursion for the moments can be readily derived using the techniques developed earlier. 



10.3 Free additive convolution of equilibrium measures 

Equilibrium measures are a fascinating topic within random matrix theory. They arise in the context of 
research that examines why very general random models for random matrices exhibit universal behavior in 
the large matrix limit. Suppose we are given a potential V{x) then we consider a sequence of Hcrmitian, 
unitarily invariant random matrices Ajy, the joint distribution of whose elements is of the form 



P{A N ) cx exp(-JVTrV(Ajv))dAjv, 
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0.7 




Figure 7: Additive convolution of equilibrium measures corresponding to potentials V\{x) and Vs(x). 

where cIAn = Y[i<j(dA.N)ij- The equilibrium measure, when it exists, is the unique probability distribution 
function that minimizes the logarithmic energy (see [11] for additional details). The resulting equilibrium 
measure depends explicitly on the potential V(x) and can be explicitly computed for some potentials. In par- 
ticular, for potentials of the form V(x) = tx 2m , the Stieltjes transform of the resulting equilibrium measure 
is an algebraic function [11, Chp. 6.7, pp. 174-175] so that the equilibrium measure is an algebraic distri- 
bution. Hence we can formally investigate the additive convolution of equilibrium measures corresponding 
to two different potentials. For V\(x) = x 2 , the equilibrium measure is the (scaled) semi-circle distribution 
encoded by the bivariate polynomial 

L mz = m 2 + 2 m z + 2. 
For V^x) = a; 4 , the equilibrium measure is encoded by the bivariate polynomial 

Lf lz = 1/4 m 2 + mz 3 + z 2 + 2/9 Vs. 

Since Ajv and B^v are unitarily invariant random matrices, if Ajy and Bjv are independent, then the limiting 
eigenvalue distribution function of C^r = Ajy+B^r can be computed from L^ z and L^ z . The limiting 
eigenvalue density function fc(%) is the free additive convolution of Ja and Jb- The Matlab command 
LmzC = AplusB(LmzA,LmzB) ; will produce the bivariate polynomial 

L„ 1Z = -9m 4 - 54m 3 z + (-108 z 2 - 36) m 2 - (72 z 3 + 72 z) m - 72 z 2 - 16 V3. 

Figure 110.31 plots the probability density function for the equilibrium measure for the potentials V\ (x) = x 2 
and V2(x) = x A as well as the free additive convolution of these measures. The interpretation of the resulting 
measuring in the context of potential theory is not clear. The matrix C^v will no longer be unitarily 
invariant so it is pointless to look for a potential V-${x) for which F c is an equilibrium measure. The tools 
and techniques developed in this article might prove useful in further explorations. 

10.4 Algebraic sample covariance matrices 

The (broader) class of algebraic Wishart sample covariance matrices for which this framework applies is 
described next. 
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Theorem 10.2. Let A n *-^A £ M. a i g , and 'Bj^i-^B G -M a i£, fee algebraic covariance matrices with G nj jy 
denoting an nx N (pure) Gaussian random matrix (see Definition \5.13\) . Let ^K. n ,N = A^G^tyB]^ 2 . XTien 

as n, N — > oo and cjv = n/iV — > c. 

Proof. Let Y n>iV = G n>N B]( 2 , T„ = Y„ iJV Y; >jv and T w = Y^Y„, N . Thus S„ = A„ x T„ = 
An T^nA-n . The matrix T„, as defined, is invariant under orthogonal/unitary transformations, though 
the matrix TV is not. Hence, by Corollary 15.271 and since A„ i— > A £ -M a ig, S„ i— > S £ M a i g whenever 

From Theorem [Ell T„nTe X alg if T w h-> f e M a i g . The matrix = B^ /2 G^ >JV G n>iV B^ 2 is 
clearly algebraic by application of Corollarv l5.27l and Theorem 15.61 since Bjv is algebraic and G^ N G n ^ is 
algebraic and unitarily invariant. □ 



The proof of Theorem 110.21 provides us with a recipe for computing the polynomials that encode the 
limiting eigenvalue distribution of S in the far more general situation where the observation vectors are 
modelled as samples of a multivariate Gaussian with spatio-temporal correlations. The limiting eigenvalue 
distribution of S depends on the limiting (algebraic) eigenvalue distributions of A and B and may be called 
using the AtimesWishtimesB function in the RMTool [22] package. See [23] for the relevant code. 

10.5 Other applications 

There is often a connection between well-known combinatorial numbers and random matrices. For exam- 
ple, the even moments of the Wigncr matrix are the famous Catalan numbers. Similarly, if Wjv(c) denotes 
the Wishart matrix with parameter c, other combinatorial correspondences can be easily established using 
the techniques developed. For instance, the limiting moments of Wjv(1) — Izv are the Riordan numbers, the 
large Schroder numbers correspond to the limiting moments of 2Wjv(0.5) while the small Schroder numbers 
are the limiting moments of 4Wjv(0.125). Combinatorial identities along the lines of those developed in [13] 
might result from these correspondences. 



11. Some open problems 

• If we are given a single realization of an N x N sized algebraic random matrix Ajy, is it possible to 
reliably infer the minimal bivariatc polynomial £ mz , i.e., with combined degree D m + D z as small as 
possible, that encodes its limiting eigenvalue distribution? 

• This is closely related to the following problem. Define the set of "admissible" real-valued coefficients 
Cij for < i < D m and < j < D z . Here admissibility implies that (a branch of) a solution m(z) 
of the equation L^ z (m,z) :— Y^ri=o CijU 1 ^ — is globally the Stieltjes transform of a positive 
probability distribution. 
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